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FORWARD 

The project goals, plan and accomplishments up to this point are summarized in the view-graphs 
following this page. Among the various accomplishments, the most important have been: 

1 . The development of the prismatic finite element code for doubly curved platforms and 
its validation with many different antenna configurations 

2. The design and fabrication of a new slot spiral antennas suitable for automobile cellular. 
GPS and PCS communications 

3. The investigation and development of various mesh truncation schemes, including the 
perfectly matched absorber and various fast integral equation methods. 

4. The introduction of a frequency domain extrapolation technique (AWE) for predicting 
broadband responses using only a few samples of the response. 

This report contains several individual reports most of which have been submitted for publication 
to refereed journals. For a report on the frequency extrapolation technique, the reader is referred to 
the UM Radiation Laboratory report 

A total of 14 papers have been published or accepted for publication with the full or partial support 
of this grant. Several more papers are in preparation. 
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LIST OF DOCUMENTS IN THIS REPORT 

A Planar Slot Spiral for Conformal Vehicle Applications 

M Nurnberger , J. Volakis, D. T. Fralick and F.B. Beck 

This document provides a brief description of the design, fabrication and measurement of the 
full scale slot spiral antenna and infinite balun feed for this antenna. The antenna covers he 
Cellular and GPS bands and is intended for automobile deployment. The measurements at the 
NASA Langley facility demonstrate the CP performance as well as improvements obtained with the 
application of a resistive coating for terminating the arms. 

Triangular Prisms for Edge-Based Vector Finite Element Analysis of Conformal 
Antennas 

T TOs papMdescribe^lK theoiy of ttie antenna analysU code FEM A-PRISM. S^vend validations 
are given todemonstrate the code's suitability for doubly curved antennas. This is the first code for 
antenna analysis on doubly curved platforms and was jointly supported by the Rome Laboratory . 

Comparison of Three FMM Techniques for Solving Hybrid FE-BI Systems 

S. Bindiganavale and J. L. Volakis , 

The paper provides a critical look at the fast multipole methods using speed and accuracy for 
benchmarking purposes. This work is towards our efforts to improve the accuracy of the finite 
element mesh truncation schemes while maintaining memory and CPU speed at tolerable levels. 

Artificial Absorbers for Truncating Finite Element Meshes 
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J Volakis, T. Senior, S. Legault, T Ozdemir and M. Casciato 

This document gives an overview of the PML absorber, its capabilities and application to 
various antenna and scattering problems, including the NASA almond and antennas on conical 
platforms 

Design of Planar Absorbing Layers for Domain Truncation in FEM Applications 

S. Legault, T. Senior and J. Volakis 

Design guidelines are given for the PML absorber. Curves and formula are derived w hich can 
be used to select the PML parameters to yield a desired performance. 

Application and Design Guidelines of the PML Absorber for Finite Element 
Simulations of Microwave Packages 

J. Gong, S. Legault, Y. Botros and J. Volakis 

Demonstrates the validity of the PML design curves derived/presented in the previous document 
for applications to microwave circuits. The latter application can be carried under a controlled 
environment and is therefore better suited for validating the PML design curves 

Hybrid Finite Element Methods for Electromagnetics: Applications to Antennas 
and Scattering 

J. Volakis, J. Gong and T. Ozdemir 

A lengthy up-to-date review of the hybrid finite element methods for scattering and 
radiation. This document is a concise look at Michigan's contributions, including mesh truncation 
schemes, feed modeling and parallelization. Many applications are included ranging from antenna 
radiation to radome performance evaluations and scattering. An introduction to the AWE frequency 
extrapolation technique is also included for modeling microwave circuits using the finite element 
method. This document was submitted for inclusion to a book prepared by ICASE. 


3 



¥ 



5 



6 



Conformal Antenna and Radome Geometries 
Simulated via the Finite Element Method 
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carried out during the 2nd year 

• New feed modeling techniques were developed in the context 
of FEM. Coaxial and aperture coupled feeds were examined 

University of Michigan 



Progress: New Archimedean Slot Spiral 

Slot Radome Modeling 
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on conformal platforms. Compared 
with measured data collected at China 
Lake(NAWC) and at MRC(Dayton) 

University of Michigan 
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coming fiscal year. 

The frequency extrapolation techniques AWE and CFH were introduced 
as a means to recover broadband responses from a few frequency data. 

University of Michigan 
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involving a mix of small and large scale details. 




A PLANAR SLOT SPIRAL FOR 
CONFORMAL VEHICLE APPLICATIONS 


M.W. Nurnberger*, J.L. Volakis*, D.T. Fralick 1 , F.B. Beck 1 


* Radiation Laboratory 
Dept, of Elect. Engin. & Comp. Sci. 
University of Michigan, Ann Arbor, MI 48109-2122 

* NASA-Langley Research Center 
Hampton, VA 23681 


ABSTRACT 

A slot spiral antenna and its associated feed are presented for conformal mounting on a variety of 
land, air, and sea vehicles. The inherent broadband behavior and good pattern coverage of the spiral 
antenna is exploited for the integration of multiple frequencies, and thus multiple transmitting and 
receiving apertures, into one compact, planar antenna. The feasibility of the broadband slot spiral 
antenna relies on the use of an equally broadband, balanced, planar, and non-intrusive feed struc- 
ture. The design of the slot spiral, its feeding structure, and the reflecting cavity are discussed with 
emphasis on the experimental validation and construction of the antenna. 

Keywords: Spiral Antennas, Planar, Conformal, Automotive, Antenna Measurements 


1.0 INTRODUCTION 

Spiral antennas are particularly known for their ability to produce very wide-band, almost perfect 
circularly-polarized radiation over their full coverage region. Because of this polarization diveristy 
and broad spatial and frequency coverage, many applications exist, ranging from military surveilance, 
ECM, and ECCM to numerous commercial and private uses, including the consolidation of multiple 
low gain communications antennas on moving vehicles. 

For the typical wire spiral antenna, the performance advantages mentioned above come at the 
price of size and complexity. While the radiating elements of a wire spiral may be planar, the feed 
network and balun structures generally are not, and combine to add weight, depth, and significant 
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compelexity to the system Furthermore, because the spiral antenna radiates bi-directiona!l\. an ab- 
sorbing cavity is generally used to eliminate the radiation in one direction, uddinc even more depth 
to the antenna. While some designs exist [ 1 ] that integrate the feed and balun into the cavitv and 
reduce the complexity somewhat, the absorbing cavity is still at least a quarter-wavelength deep ai 
the lowest frequency of operation, adding significant thickness to the antenna [2], 

A slot spiral is not burdened with most of these difficulties. As was demonstrated previously [3], 
the balun and feed structure may be integrated into the planar radiating structure. This greatly sim- 
plifies the construction and increases the accuracy of the antenna by allowing standard prmted circuit 
techniques to be used throughout the entire facrications process. Also, while the slot spiral also radi- 
ates bi-directionally, a deep, absorbing cavity is not necessary for uni-directional radiation. In fact 

for conformal mounting, a very shallow reflecting cavity is sufficient, and also serves to increase the 
gain somewhat. 


2.0 SPIRAL ANTENNA OPERATION 

The operation of a spiral antenna may most easily be understood by considering a two-wire transmis- 
sion line that has been wrapped around itself to form a spiral. Clearly, prior to its deformation, the 
two-wire line did not radiate, as the currents in the adjacent wire were out of phase, each cancelling 
the other s radiation in the far field. Following the rearrangement, at the center of the spiral the 
currents on the two adjacent wires are still 180° out of phase. However, as the currents travel down 
the two wires, because they have traveled different distances, they reach a region where instead of 
being out of phase with their neighbors, they are in phase. This annular region is called the active 
or radiating region, because here, instead of cancelling the effects of their neighbors, each current 
reinforces the others, creating appreciable radiation in the far field. It is interesting to note that, for 
each current element in this active region, there is another current element that is in both space and 
phase quadrature. Thus, the radiation from all the current elements in the active region combines in 
the far field to produce circularly polarized radiation. 

In order to be considered frequency independent, an antenna must obey both the “angle princi- 
ple” and the “truncation principle,” as defined by Rumsey [4], Unfortunately, from its geometry, and 
as can be concluded by studying Figure 3, the Archimedean spiral antenna follows neither Figure 
3 shows the radiation pattern of a spiral antenna, operated at a frequency inside its “frequency inde- 
pendent” frequency range, measured with a spinning linear source. Specifically, the high axial ratio 
indicates that there is a significant amount of opposite sense circularly polarized energy combining 
with the desired sense circularly polarized radiation. This behavior, which is present over the entire 
operational frequency range, has been determined to be a result of the truncation of the spiral arms 
[2,5], and indicates that the currents are not dying out before reaching the end of the antenna. Rather 
they are traveling to the end of the antenna, and reflecting back along the spiral arms into the active 
region, from which they re-radiate energy with the opposite sense circular polarization. Thus, a large 
portion of the effort in designing an Archimedean spiral antenna lies in minimizing this reflection. 
In essence, the structure must be made “frequency independent.” 

Another important aspect of spiral antenna design is the necessity of a “frequency independent ” 
balanced feed. Several different techniques have been developed [1,2], most notably Dyson’s “infi- 
nite balun [6], a form of which is used in the antenna under development. Besides providing a bal- 
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anccd feed from the unbalanced coaxial line, a properly designed balun prevents the pattern squint 
so often observed in early spiral antenna designs. 

Most spiral antenna designs also include a backing eas ily to force uni-directional operation. It 
has been determined that, for wire spirals, an absorbing eas ily offers the best trade-off in terms of 
frequency bandsvidth and physical size [2.5], For slot spirals, hosvever. a reflecting cavity can be 
used to great effect, as the available bandsvidth is much larger than that obtained when coupled svith 
a wire spiral. 

From the abos'e description it seems reasonable, and in the follosving discussion it will prose 
invaluable, to think of the spiral antenna in terms of three different regions of operation — the trans- 
mission line region, svhere very little radiation takes place, the active region, svhere the predominate 
amount of energy is radiated, and the termination region, where the unradiated energy is absorbed. 
This division of the antenna into smaller, almost independent pieces allows certain parameters and 
effects to be separated out and better understood by themselves before trying to understand the op- 
eration of the antenna as a whole. 

The above discussion has. for clarity’s sake, implicitly assumed a wire spiral, and its associated 
electric currents. It may be extended to include the slot spiral by applying the theory of duality, in 
which the two-wire line is replaced with a slot line, and the electric currents with magnetic currents. 


3.0 DESIGN INFORMATION 

The design procedure for the slot spiral antenna assembly is somewhat long and involved, and re- 
quires the careful choice of many inter-related design parameters. In this section we summarize the 
concerns in the design of the slot spiral, and discuss the more important design parameters briefly. 
Some guidelines are also given to assist the designer in developing a useful antenna. 

The design parameters are: 

• spiral growth rate, a 

• number of turns, N 

• microstrip line characteristic impedance(s), Z M . 

• slot line characteristic impedance(s), Z s , 

• substrate characteristics: e r , tan 6, metalization technique, ... 

• thickness of the substrate, t 

• length and impedance of the tuning stub, l Stub & Z Stub 

• cavity depth, h 

• cavity wall material 

• slot line termination material and technique 
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• specific geometry of the microstrip-slot line transition 

As with most antennas, the design procedure of the .slot spiral is an iterative one. However, be- 
cause the individual parameters are highly dependent on each other, the design procedure is particu- 
larly involved. It is perhaps easiest initially to choose the growth rate, and then to carefully inspect 
the initial design for inconsistencies. For the following discussion, please refer to Figures 1 ** and 
4-6 for clarification. 

In the 500 MHz to 4 GHz region, the authors have found a growth rate of 65.4 mils/rad to be a 
good initial choice. To provide a good match at the connector, the initial characteristic impedance 
of the microstrip line is chosen to be 50.Q. The characteristic impedance of the m.crostr.p inside the 
spiral will be decided later, based on the inter-slot spacing and construction feasibility. At the center 
of the antenna, the actual feed is contrived through the use of a microstrip to slot 'line transition 
Here, the microstnp line views the slotline as a pair of shunt branches. For the transition to operate 
most effectively, the characteristic impedance of the microstrip line is chosen to be half that of the 
slot line, yielding a perfect match at the feed. Based on the dielectric constant and thickness of the 
substrate, and on fabrication concerns, a certain range of slot line impedances will be available For 
a dielectric constant of 4.5 and a thickness of 15 mils, the authors have found that a 100fi slot line, 
which has a width of approximately 20 mils, is convenient. This yields a microstrip characteristic 
impedance of 50fi, and a microstrip groundplane width of approximately 185 mils. 

To minimize coupling between the mictrostrip feed line and the slotline, it is best to maximize 
the ratio of microstrip line to ground plane width. Many design issues can be traded off in this case 
— the most important are the growth rate, the substrate thickness and dielectric constant, and the 
slotline and microstnp characteristic impedances. In this iteration just the microstrip line width will 
be modified — the others may be adjusted in later iterations. To maximize the width ratio mentioned 
above, the impedance of the microstrip line is increased until fabrication difficulties determine the 
minimum width. For safety, the authors don’t venture below a 5 mil feature size. Klopfenstein ta- 
pers [7] are used in both cases to ensure that these impedance changes do not hinder the broadband 
behavior of the antenna and feed. 

To terminate the microstrip line following the transition, for expediency the authors currently use 
a Xg/A open-circuit stub of significantly lower impedance than the microstrip line. More wide-band 
techniques are currently being investigated. The termination of the slot line, as discussed above is 
critical for making the antenna behave in a frequency independent manner, and is the subject of on- 
going investigation. Initially, tapered applications of foam and ferrite-loaded rubber absorbers were 
used, but were found to be of little help. Currently under investigation is the use of “radar-absorbing 
paint, applied in a tapered fashion over the slots. 

The design of the reflecting cavity is driven by the upper and lower frequency extents of the 
antenna. The reflector must be close enough to the antenna so that, at the highest frequency, it is still 
less than A 0 /4 away, to avoid complete cancellation on boresight. However, it cannot be so close 
that it shorts the fields in the transmission line region. In practice, these fields are very tightly bound 
to the slots, but it is still quite possible to significantly disturb them. To help avoid resonant cavity 
effects, the walls of the reflecting cavity are also made of absorber. 

As was mentioned earlier, the quality of the slot line termination is very important for low axial 
ratio designs and consistent pattern quality over frequency. The growth rate is also important, and 
should be kept as small as possible, again to minimize the axial ratio. The microstrip to slot line 
transition is important as well, particularly in terms of the overall bandwidth and impedance match of 
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the antenna. Finally, choice of the substrate is critical to minimize losses, especially in the microstrip 
line, and to reduce unwanted coupling effects 


4.0 MEASUREMENTS 


A slot spiral antenna, with parameters as discussed above, and operating frequency range from 750 
MHz to 1.8 GHz, was conducted on 15 mil Rogers TMM-4, with 1 oz. electrodeposited copper on 
each side. Figures 4-6 show the front, back, and disassembed views of the antenna itself. Calibrated 
pattern measurements were conducted at the NASA LaRC low frequency antenna chamber, using a 
spinning linear source antenna. The antenna was mounted in a 4.5’ square ground plane attached to 
the roll-over- azimuth-over-elevation positioner in the chamber, as depicted in Figure 7. The edges 
of the ground plane were also treated with tapered R-card to reduce diffraction effects. 

Orthogonal plane cuts were taken through and at 90° to the connector in an effort to determine 
and minimize radiation effects from the connector. Figure 8 show's a typical pattern cut, orthogonal 
to the plane of the connector, with a plot of the axial ratio as well. This pattern was taken with the 
slot line termination constructed from ferrite-loaded rubber and absorbing foam. Figure 9 shows the 
radiation pattern for the same antenna, but with a tapered application of “radar-absorbing paint” as 
the termination. The decrease in axial ratio is significant, and demonstrates the necessity of a good 
termination. Patterns similar to these were obtained over the entire operating range of the antenna, 
each showing a similar reduction in axial ratio. 

Figures 3, 8, and 9 also show an unfortunate lack of gain — the peak gain is approximately -2.5 
dBic. Preliminary calculations indicate that most of this loss is attributable to losses in the microstrip 
line, especially since the copper on the substrate was electrodeposited. Techniques to minimize these 
losses are currently being studied. 


5.0 CONCLUSIONS 

The design of the slot spiral antenna was discussed, and experimental results were obtained to exam- 
ine its performance. A useful broadband slot termination is now being investigated with promising 
results, implying that the Archimedean slot spiral can be made truly frequency independent. Cur- 
rent and future work includes further numerical and experimental studies of the slot line termination, 
synthesis of a more broadband microstrip to slot line transition, and pattern performance evaluation 
when mounted on automobiles or other moving vehicles. 
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Spiral Sloiline Spiral Microstrip Balun/Fced 


Figure 1 : Schematic diagram of the slot 
spiral antenna and microstrip 
balun/feed. 



Figure 2: Enlarged view of the Section A- 
A\ Figure 1, showing the antenna 
cross-section and one slot termina- 
tion technique. 
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Azimuth (deg) 


Figure 3: Spinning linear plot of spiral antenna, showing effects of poor arm termination 
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Figure 4: Top view of Slot Spiral Antenna, 
showing infinite balun microstrip 



Figure 5: Bottom view of Slot Spiral An- 
tenna, showing absorber walls and 
reflecting cavity bottom. 
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Figure 6: Disassembled pieces of the Slot 
Spiral Antenna. 


Figure 7: Slot Spiral Antenna mounted on 
ground plane, in anechoic chamber. 
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Figure 8: Spinning linear plot of spiral antenna, with ferrite-loaded rubber and absorbing foam 
termination. 
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Figure 9: Spinning linear plot of spiral antenna, with “radar-absorbing paint” arm termination. 
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Abstract 

This paper deals with the derivation of the edge-based shape functions for the 
distorted triangular prism and their applications for the analysis of doubly curved 
conformal antennas in the context of the finite element method (FEM). Although 
the tetrahedron is often the element of choice for volume tessellation, mesh gen- 
eration using tetrahedra is cumbersome and CPU intensive. On the other hand, 
the distorted triangular prism allows for meshes which are unstructured in two 
dimensions and structured in the third dimension. This leads to substantial sim- 
plifications in the meshing algorithm and many conformal printed antenna and 
microwave circuit geometries can be easily tessellated using such a mesh. The 
new edge-based shape functions are first validated by computing the eigenvalues 
of three different cavities (rectangular, cylindrical and pie-shell). We then proceed 
with their application to computing the input impedance of conformal patch anten- 
nas on planar, spherical, conical and other doubly curved (ogival) platforms, where 
the FEM mesh is terminated using an artifical absorber applied conformal to the 
platform. Use of artificial absorbers for mesh termination avoids introduction of 
Green’s functions and, in contrast to absorbing boundary conditions, a knowledge 
of the principal radii of curvature of the closure’s boundary is not required. 

^his work was supported in part by the U.S. Air Force Rome Laboratory and the 
NASA Langley Research Center. 
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1 Introduction 


The brick and tetrahedron are popular elements for a finite element analy- 
sis of electromagnetic problems. The first is indeed attractive because of its 
simplicity in constructing volume meshes whereas the tetrahedron is a highly 
adaptable, fail-safe element. It is often the element of choice for three dimen- 
sional (3D) meshing but requires sophisticated and CPU-intensive meshing 
packages. The distorted prism (see Figure 1) is another volume element 
which provides a compromise between the adaptability of the tetrahedron 
and the simplicity of the brick. Basically, the distorted prism allows for 
unstructured meshing (free-meshing) on a surface and structured meshing in 
the third dimension. An approach for growing prismatic meshes is illustrated 
in Figure 2 and most volumetric regions in antenna and microwave circuit 
analysis can be readily tessellated using such a mesh. As seen, once the sur- 
face grid over the platform is constructed, the volume mesh is grown along 
the surface normal by repeating the same grid at multiple distances from 
the platform. This avoids use of CPU-intensive volumetric meshing packages 
and, in many cases, including some popular patch shapes, the surface grid 
can be easily generated without resorting to a surface gridding package. Ex- 
amples include printed rectangular and circular patch antennas, and circuits 
comprised of rectangular shapes. Moreover, because of their triangular cross- 
section, the prisms overcome modeling difficulties associated with bricks at 
corners formed by planes or edges intersecting at small angles. 

A special case of the distorted prism is the right prism which is character- 
ized by the right angles formed between the vertical arms and the triangular 
faces [1]. The top and bottom faces of the right prism are necessarily parallel 
and equal, restricting them to a limited range of applications, namely, ge- 
ometries with planar surfaces. In contrast, the distorted triangular prism is 
almost as adaptable as the tetrahedron with the exception of cone-tips which 
are not likely to occur in printed antenna and microwave circuit configura- 
tions. 

In this paper, we introduce edge- based shape functions [2]- [4] for the most 
general distorted prism shown in Figure 1. These prisms have non-parallel 
triangular faces and each of their three vertical edges can be arbitrarily ori- 
ented. In the following, we first present the derivation and validation of 
the edge-based shape functions and then proceed with their application for a 
characterization of printed antennas on doubly curved surfaces. In this appli- 
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Figure 1: The distorted triangular prism shown with the direstions of the 
edge vectors 

cation, the finite element mesh is terminated conformal to the platform's sur- 
face using an artificial absorber rather than an absorbing boundary condition 
(ABC) [5], The employed conformal mesh termination is easily implemented 
by using prisms, and in contrast to the ABCs, the artificial absorber does not 
require a prion knowledge of the closure’s radii of curvature or the wave's 
propagation characteristics. The utility and versatility of the proposed fi- 
nite element method (FEM) formulation is demonstrated by considering the 
analysis of several printed antennas on different platforms. Specifically, we 
include input impedance computations for rectangular and circular patches 
on planar, spherical and conical surfaces. The radiation from a patch on an 
ogive is also considered. 

2 Vector edge-based basis functions 

Consider the distorted prism shown in Figure 3. The prism’s top and bot- 
tom triangular faces are not necessarily parallel to each other and the three 
vertical arms are not perpendicular to the triangular faces. The first step 
toward specifying a set of shape functions for the prism is the identification 
of a unique cross-section containing the observation point (x,y,z) (see Fig- 
ure 3). This is done by introducing a parametric representation for the nodes 
(x,, y u z t ) of this cross-section as illustrated in Figure 4. These parametric 
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Figure 2. Prismatic mesh: Unstructured over the antenna surface, structured 
along the surface normal. 


Perspective view 


Top view 



> \ 

---> N Interior point (x,y,z) ’ 


Triangular cross-section 
containing the interior point 

Figure 3: Variation of the triangular cross-section along the length of the 
prism 
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Figure 4: (a) Nodal coordinates, (b) triangular cross-section with the local 
coordinates £ and rj. 


representations are given in the Appendix. They involve the parameter 5 
such that (x,,y,,Zi) reduce to the nodes of the bottom triangle when s = 0 
and those of the top triangle when 5 = 1. Given a point (x,j/,r) interior to 
the prism, a unique value for s can be computed as discussed in Appendix A. 

Having specified the cross-section through the point (x,y,z), we next 
proceed with the derivation of the basis functions. We chose to represent 
the field variation across the triangular cross-section using the Whitney- 1 
form [6]. A simple linear variation will be assumed along the length of the 
prism. Specifically, the vector basis functions for the top triangle edges can 
be expressed as 


Nl = d\ ( L* 2 ^ — Z13VZ/2) 5 

N2 = d .2 ( Z/3VZ1 — L1VT3 ) s (1) 

N3 = <f 3 ( L1VL2 — L 2 ^Li\)s 

and correspondingly those for the bottom triangle edges will be 

= d.\ ( Z/2 V Lj — T3 V L2 ) (1 — s) 

M 2 = <^2 ( T 3 V L\ — L\ VT 3 ) (1 — s) ( 2 ) 

M 3 = <^3 ( L\ V Z *2 — L 2 V L \ ) (1 — s). 

The subscripts in these expressions identify the edge numbers as shown in 

Figure 1 and the distance parameters d, are equal to the side lengths of the 
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(a) 



Figure 5 : (a) \ector map of N 2 or M 2 , (b) Vector map of (c) Variation 
ot v as a function of x , y, and z. 


triangular cross-section containing the observation point (see Figure 4(b)). 


L\{Z,V) 

. l 3 U,n) 

are the usual two dimensional scalar node-based basis functions [7] for the 
same triangle with a, denoting the interior vertex angles and h t being the 
node heights from the opposite side. The variables { and r) represent the 
local coordinates and are illustrated in Figure 4(b). As required with all 
edge-based shape functions, N, • e, and M; • e; have unity amplitude on the 
Ith ed g e whereas N, • e, = M, • e, = 0 for * ^ j. Their vector character is 
depicted in Figure 5 (a) and they simply ”curl” around the node opposite to 
the edge on which their tangential components become unity. 

It remains to define the shape functions for the three vertical edges and we 


cos 03 sina 3 
cosa 2 sino 2 

~r-Z + — L — »?• 


( 3 ) 
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chose to express these by the representations (linear over the cross-section i 


Ki( s f .r,) = v(‘.r))L x (i. tj) 

t ?) = i'it.ri) L 2 U.t]) ( 4 ) 

K 3 (^, 77) = r(£,7?) LsiZ.T}). 


As before, L x are the node-based shape functions defined in (3) and a pictorial 
description of K! is found in Figure 5(b). Of particular importance in (4) is 
the unit vector v(£.r}). It is a linear weighting of the unit vectors tq, v 2 and 
r 3 associated with the vertical arms (see Figure 5(c)). and is given by 


$(£.»?) 


E?=l LitfitfVj 

11 uLi *?)£’. 1 r 


(5) 


This particular choice of v is oriented parallel to the side faces of the prism 
when evaluated on those surfaces and minimizes tangential field discontinuity 
across the side faces. Another choice is 


v{(,v) = Vs ( 6 ) 

which is always oriented normal to the triangular cross-sections of the prism 
and ensures tangential field continuity across the top and bottom triangular 
faces. Both choices are equally useful. However, (5) is more computationally 
efficient and has been used in generating the results presented in sections 3 
and 4. 

The shape functions derived above for the distorted prism reduce to the 
right prism shape functions presented by Nedelec [8]. However, it should be 
pointed out that the above shape functions do not guarantee tangential field 
continuity across the faces of neighboring prisms. The possible discontinuity 
is primarily an issue for highly distorted elements and is caused by the fact 
that the horizontal vector basis functions (N and M) may have small non- 
vanishing tangential conponents across the vertical faces of the prism (and 
vice versa for the vertical basis functions). The expressions given in Ap- 
pendix B neglect contributions due to the tangential discontinuities across 
the inter-element boundaries. These contributions were ignored because the 
discontinuity depends on the distortion of the prism and in practice the ele- 
ments are marginally distorted. For general applications, our basis functions 
can be safely employed (within the accuracy of the finite element method) as 
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long as the deflection of the vertical arms is less than o' 5 . In the results to be 
presented later, the deflection was always kept less than 3'. Furthermore, in 
the case of patch antennas, the discontinuity is nearly nonexistant because 
the horizontal prism faces (top and bottom) coincide with metallic surfaces 
and therefore any possible discontinuity on the vertical fields under the patch 
(which are most dominant) is eliminated. 

Shape functions ensuring tangential continuity across the distorted prism 
faces can be derived but will be more complicated in form and CPU intensive 
to implement. On the contrary, the ones presented here lead to a more useful 
code for design purposes. 

3 Eigenvalue Computation 

In this section, we examine the validity of the presented edge-based functions. 
Specifically, we consider the eigenvalues of three different cavities using the 
distorted prism as the tessellation element. We begin by first deriving the 
matrix elements following Galerkin’s testing. The weighted residuals of the 
vector wave equation are 

/ I Jv N ' * V X E ~ k o E ) dv = 0, i = 1,2,3, (7) 

/ / ^ M t - • (V x V x E - ^ E)dV = 0, * = 1,2,3, (8) 

/ //K,-(VxVxE-*’E)dV = °, i = 1,2,3, (9) 

in which N,, M, and K, comprise the nine edge- based vector basis functions 
defined in the previous section and E is the electric field vector. 

The matrix equations are generated by introducing the representation 

3 

E = ^[Eis N,-(r) + EiMir) + £^K t (r)] (10) 

»=i 

where Eis , EiM and Eix are the expansion coefficients, and correspond to 
the average amplitudes of the tangential field vector along the prism’s edges. 

Substituting (10) into (7)-(9), and invoking the divergence theorem yields 
the element equations 

E EjslNNCij - k]NNDij\ + E E ]M [NMC t] - k]NMD tJ ] + 

]= l j = 1 
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3 

E Ej K [ .V A C,J - k; X A D, , ] = 0 i 1 1 1 

J = 1 

E Ejx[MXC tJ - klMXD tJ ] + £ Ejm[M MC, j - k 2 0 MMD t] } + 

;=i J = 1 

- k 2 0 MKD t] ] = 0 (12) 

;= i 

£ £>[A\VC 0 - klKND,,} + £ £ jJf |A-MC„ - tjAitfA,] + 

J=1 J=1 

3 

X] AjA [A' A' C 0 - k 2 0 KI\ D tJ ] = 0, t = 1, 2, 3. ( 13) 

j=i 

where the quantities NNC,j, NND etc. are integrals (over the volume 
of the prism) involving vector basis functions. Explicit expressions for the 
integrals are given in Appendix B. 

Upon assembly of (11)- (13) and boundary condition enforcement, we ob- 
tain the generalized eigenvalue system 

[A]{x} = k][B){x] (14) 

in which A = k 2 represent the eigenvalues of the problem. The matrices [A] 
and [B\ are real, symmetric and sparse ([j 9] is also positive definite). 
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Figure 6. (a) Rectangular cavity discretized using right triangular prisms, 
(b) eigenvalues for the air-filled cavity and comparison with the brick and 
tetrahedron results. 

Example 1: Rectangular Cavity 

The first example is the rectangular cavity shown in Figure 6(a). Results 
based on brick, tetrahedron [9] and prism discretizations are given in Fig- 
ure 6(b), where they are compared to the exact values. We observe that, 
overall, the data based on the prism are better than those based on the brick 
and, at least, as good as those associated with the tetrahedron. Figure 6(a) 
displays the actual mesh used for both the brick and prism discretizations. 
The number of degrees of freedom for the prism, brick and tetrahedron com- 
putations were 382, 270 and 260, respectively. 

Example 2: Dielectric ring resonator 

Figure 7a shows a dielectric ring placed symmetrically inside a cylindrical 
cavity and of interest is the computation of the resonance frequency in the 
presence of the dielectric ring. The resonance frequency was measured [10] 
using a loop antenna connected to a directional coupler as shown and the 
cavity was excited by the same loop antenna. Since maximum coupling to 
the cavity occurs at resonance, minimum power is returned to the detector 
at this frequency. For computation, the cavity was discretized as shown in 
Figure 7b resulting in 1051 degrees of freedom. The computed frequency was 
1257 MHz (based on the smallest eigenvalue) and this is within 2% of the 
measured value (1282 MHz). 
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Dielectric ring 



Cavity dimensions : 
Height = 2.2cm 
Inner diameter = 7.5cm 

Outer diameter = 1 1cm 
E r = 1 3.5 


Dielectric ring info : 
Height = 13.84cm 
Diameter = 15.24 cm 


(b) 

Results. 

Measured = 1282 MHz 
Computed = 1257 MHz 
( Error = 2 %) 


Figure ?: Dielectric ring resonator, (a) geometry, (b) mesh. 


Example 3: Pie-shell Cavity 

The third example is a pie-shell sector as shown in Figure 8(a) modeled 
using distorted prisms. It is obtained by bending the rectangular cavity 
considered earlier. The computed and exact eigenvalues for the first five 
dominant modes are given in Figure 8(b), and these testify to the accuracy 
of the distorted prisms in modeling curved geometries. The number of degrees 
of freedom used for this computation was 382. 

4 Application to antennas 

4.1 Finite element-artificial absorber method 

Prisms facilitate the modeling of conformal antennas since the presence of 
curvature does not complicate the mesh generation. As is the case with all 
PDE methods, the mesh must be terminated using a boundary integral, some 
approximate boundary condition or an absorber. When a Green’s function is 
available, the mesh termination can be achieved using the boundary integral 
method right on the antenna surface yielding exceptionally accurate results. 
However, there are two problems with the boundary integral mesh trunca- 
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Figure 8: (a) Pie-shell cavity discretized using distorted prisms, (b) eigen- 
values for the air-filled cavity. 
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Figure 9: Finite element-artificial absorber technique: (a) Antenna geometry, 
(b) Mesh termination using a metal-backed lossy dielectric layer. 
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tion techniques. First, the Green's function is only available for canonical 
geometries, i.e.. planar [11], cylindrical [12] and spherical [13]. Second, the 
presence of material overlays can prohibitively complicate the derivation and 
evaluation of the Green’s function making it very difficult to implement the 
boundary integral termination method. Thus, for modeling non-canonical 
geometries with possible material overlays, it is preferable to avoid use of 
the Green s function. Instead, an artificial absorber will be used for truncat- 
ing the mesh as illustrated in Figure 9. The resulting implementation will 
be referred to as the finite element-artificial absorber (FE-AA) method and 
has several attractive computational features. Among them are fast conver- 
gence [5] and a capability to provide accurate input impedance computations. 
Accurate radiation pattern data can also be obtained when the rest of the 
platform has little or no effect. In its simplest form, this is done by ignor- 
ing the extent of the platform and integrating the aperture fields using the 
free-space Green's function. 

We next discuss the artificial absorber and then proceed with the appli- 
cation of the FE-AA method for the analysis of a variety of patch antennas 
on different platforms, some of which are presented here for the first time. 

4.2 Design of the artificial absorber 

The absorber consists of a lossy dielectric layer backed by a metal. Metal 
backing enables us to terminate the mesh while the lossy dielectric lining 
traps the incoming wave and absorbs it, thereby forming a trasparent bound- 
ary. The absorber material parameters are chosen to completely eliminate 
wave reflections at normal incidence. Away from normal, the absorber reflec- 
tivity increases monotonically but can be minimized by proper selection of 
the absorber thickness and material parameters. Clearly, a thicker absorber 
has a better performance but carries additional computational burden. In 
addition, higher attenuation is achieved by making the layer more lossy. 
However, in this case more samples are required along the thickness of the 
layer to accurately model the field’s amplitude decay. For a given thickness 
and sampling, an optimization can however be carried out. Such a study was 
performed in [5] and it showed that a minimization of the reflectivity for a 
3-layer, 0.15A o (free-space) thick metal-backed absorber yields the constitu- 
tive parameters of e r = = 1 — j 2.7 (see Figure 9). Below, we make use 

of this absorbing layer for mesh truncation. A layer based on the recently 
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Figure 10. Cavity-backed circular patch on planar platform, (a) Configu- 
ration, (b) Input impedance vs frequency and comparison with the finite 
element- boundary integral method. 


introduced anisotropic perfectly matched absorber [14] could have been used 
as well, but this was not found necessary at this stage. 

4.3 Input impedance computations 

Cavity-backed, circular patch on planar platform 

We first consider a patch on a planar platform since validation data can be 
generated using the more rigorous finite element-boundary integral method. 
Figure 10a shows a circular patch backed by a circular cavity recessed in a 
planar ground plane. A probe feed has been placed between the patch and 
the ground plane and Figure 10b shows the variation of the input impedance 
as a function of the excitation frequency. As seen, the computed impedance 
is in excellent agreement with the reference data obtained by terminating 
the finite element mesh at the cavity aperture using the half-space Green’s 
function [11]. 
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For this example, the computational domain was discretized using 7.440 
prisms resulting in 12.523 degrees of freedom. One frequency run took 3.5 
minutes on an HP 715/64 workstation. 

Microstrip circular patch on sphere 

Figure 11a shows a microstrip circular patch placed on a sphere. This is 
an example which clearly shows the advantage of the finite element-artificial 
absorber technique. Once the antenna aperture is triangulated, the volume 
mesh is grown along surface normals (using prisms) and terminated with 
the artificial absorber. Figure lib shows the comparison with the moment 
method [13] and the effect of sphere radius on the resonance behavior. Res- 
onance frequencies predicted by this method and the moment method are 
within 2% of each other. Although negligible, the difference must be evalu- 
ated in view of the fact that the reference data is based on the approximation 
that only the dominant T M n mode exists under the patch. Such assump- 
tions are absent in the finite element analysis. Clearly, the value of the input 
resistance is a strong function of the employed feed model and therefore, it 
is not surprising to see differences in the levels of the resistance as computed 
by the finite element method and the single mode moment method solution. 
Figure lib also shows the resonance behavior for different sphere radii, and 
it is observed that the resonance frequency decreases with increasing radius. 
This trend might be explained by noting that the shortest distance between 
the radiating edges of the patch is greater for a larger sphere radius. The 
patch radius (measured on the spherical surface) has been kept constant for 
this calculation. 

Figures 12c shows the input resistance variation with the excitation fre- 
quency for different substrate/superstrate material constants and thicknesses. 
Similarly to the planar patch, we observe the downward shift in the resonance 
frequency caused by the increasing relative permittivity of the substrate. No- 
tice also that the shift is only half as much when the superstrate is present 
even though the increase in the refractive index of the superstrate layer is 
1.6 times higher than that of the substrate. This is because the field is much 
stronger under the patch than above it. As expected, Figure 11c shows that 
the loss in the substrate material has no effect on the resonance frequency 
but it lowers the overall level of the input impedance. The typical bandwidth 
increase with increasing substrate thickness is also observed. 
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Figure 11: Microstrip spherical-circular patch, (a) Geometry, (b) Comparison 
with reference data and effect of sphere radius on the resonance behavior, 
(c) resonance behavior as a function of substrate/superstrate thickness and 
material constants, (d) definition of the parameters. 
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Figure 12: Microstrip sectoral patch on cone: (a) Configuration, (b) change 
in the resonance behavior as a function of the cone angle 6. 


For this example, the computation region has been discretized using 6, 048 
prisms resulting in 9,997 degrees of freedom. One frequency run took 4.5 
minutes on an HP 715/64 workstation. 

Microstrip sectoral patch on cone 

This is an example which illustrates the effectivenenss of the new tech- 
nique to model antennas on doubly curved platforms with varying radii of 
curvature. To our knowledge, this is the first analysis of patches on such a 
platform. Figure 12a displays a sectoral microstrip patch printed on a coni- 
cal ground plane and Figure 12b shows the antenna resonance behavior as a 
function of the cone angle 6. We clearly observe that the resonance frequency 
drops with the cone angle for the same patch dimensions. It should be also 
remarked that the computed resonance frequency is within 3.2 percent of 
that predicted by the approximate cavity model, and this is reasonable. 

In generating the results given in Figure 12, the computational domain 
was discretized using 2,358 prisms resulting in 3,790 degrees of freedom. 
One frequency run took 5.5 minutes on an HP 715/64 workstation. 
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Figure 13: Cavity-backed rectangular patch on ogive: (a) Setup, (b) antenna 
dimensions, (c) Comparison with the measurement. 


4.4 Radiation pattern computations 

Cavity-backed rectangular patch on ogive 

As another example, we present the FE-AA method’s application to a 
more general doubly curved platform. Figure 13a shows the set-up, where a 
rectangular patch is placed on the aperture of a rectangular cavity recessed 
in the ogive’s surface. Also, in Figure 13b we show the dimensions of the 
antenna. The radiated field is computed by integrating the tangential aper- 
ture fields using the free-space Green’s function. Thus, the shown pattern 
accounts only for the antenna (and the curvature of the platform), but does 
not include interactions with the ogive’s tips. Figure 13c shows the com- 
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puted ^-polarized radiation as compared to the measurement [15]. Clearly, 
the agreement is very good for this polarization. However, predicting the o- 
polarized radiation (not shown) requires modeling of the entire ogive as this 
particular polarization has a vertical surface field component which is known 
to cause diffraction from the ogive's tips. A way to account for such sec- 
ondary diffractions is to interface the FE-AA method with a high frequency 
technique and an encouraging study in this direction has been carried out in 
[16]- 

5 Conclusions 

A new finite element technique was introduced for the analysis of printed 
antennas recessed in platforms of non-canonical shape. The distorted prism 
was introduced as the volume discretization element, and corresponding edge- 
based shape functions were derived and tested for eigenvalue computations. 

A major part of the paper though was devoted to characterizations of 
printed antennas using the new technique. Use of prismatic elements was 
found very attractive for this application and was essential in simplifying the 
modeling of antennas on doubly curved platforms. An artificial absorber was 
used to terminate the mesh conformal to the platform thereby minimizing the 
size of the problem while preserving the sparsity of the finite element matrix. 
Use of the absorber also avoids difficulties associated with the conformal 
application of the classical ABCs. Antennas on spherical, conical and ogival 
platforms were considered without using a Green’s function and therefore, 
the superstrate/substrate metarials can be readily accounted for. 

A limitation of the technique is that it models the immediate neighbor- 
hood of the antenna, and therefore ignores the details associated with the 
rest of the platform and possible substructures around the radiator. How- 
ever, the proposed method was shown to be a good approach for predicting 
the resonance behavior of antennas, and could evolve as an important tool 
for designing conformal antennas on doubly curved platforms. 
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Appendix A 

Referring to Figure 4a. a way to specify the nodes of the triangular cross- 
section through the point (x.y.z) is to use the parametric representation 

r : = r,<, + s(r lt — r.f,). i= 1.2.3. (.41) 

" here x, are the nodal position vectors of the triangular cross-section (see 
Figure 4a) whereas r„ and r lfj are associated with the top and bottom tri- 
angular faces, respectively, and 0 < s < 1. Thus, for s = 0 r, specify the 
nodes of the bottom triahgle and for s = 1, r, reduce to the nodes of the 
top triangle. The cross-section sweeps the entire volume of the prism as s 
assumes values between these two limits. 

Given an interior point (x, y , z), the task is to specify the cross-section of 
the prism that contains the point. In mathematical terms, a solution of s in 
terms of x,y, and z is required. On our way to finding this solution, we note 
that the normal distance of the point (x,y,z) to the cross-section must be 
zero, viz. 


{x - Xi)n x + (y - y x )n y + (z - z x )n z = 0 (.42) 

where x u y x and z x are the coordinates for one of the nodes of the cross- 
section (see Figure 4a), and n r , n y and n z are the components of the cross- 

section s unit normal pointing toward the top surface and can be computed 
as 


nfs) = ( r 2~ r i) x ( r 3-f2) 

II (r 2 — r x ) X (r 3 — r 2 ) || 

where r, are as given in (Al). The left hand side of (A2) is a linear function 
of s for the right prism, and is nearly linear for a marginally distorted prism 
(which is the case in practice). Therefore, a very fast converging solution of 
(A2) is obtained by the Newton- Raphson method. In generating the results 
presented in sections 3 and 4, an average of three iterations sufficed to achieve 
a tolerance of 10 -5 x (Phsm Height). 
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Appendix B 


For the assembly of the finite element equations, the following quantities 
must be computed: 


\xc t( 
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where the integrals axe over the volume of the prism and they must be eval- 
uated numerically. However, numerical integration over the distorted prism 
volume is cumbersome and therefore, the distorted prism is first mapped onto 
a right prism as shown in Figure Bl. The integration over the new volume is 
then carried out by using the five point Gaussian formula. The relationship 
between the two integrations is given by 


I f(x,y,z)dxdydz = / / / f(x,y,z)det[J]d(drjd^ 

J J JV X y Z J J 
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where (x.y.~) are related to (s -7-0 through the bilinear transformation 
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and (x,6, 3/»* ) and (x, t , y tt , z it ) are the coordinates of the prism’s bottom 
and top surface nodes, respectively, as shown in Figure Bl. Also, [J] is the 
Jacobian of the bilinear transformation defined by (Bl)-(B3) and is given by 


[J} = 


dx 

djt 

dz 

f« 

di 

d£ 

dx 

ski 

dz 
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dr) 

dx 


dz 

d< 

dC 
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Special Case: Right Prism 


For the right prism, closed form expressions can be derived for the above 
integrals. Referring to Figure 4b, we have 


ENNCit 


djdt , cos 0kn cos 8 jm cos fi km 

c k,K X,m hjh m 

cos $jn 2 c* h\d\ 

Tk xtm + 3 wX sinSik 3,n > > 


Xjn 


46 




(Distorted prism) (Right prism) 


Figure B 1 : Mapping distorted prism onto a nght prism 
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E.XKD,/ = 0 
EMM D tl = EXXDu 
EM I< Du = 0 


EKKDu = c\u 
where c is the height of the prism, 

$ rs - I 0 if r = 5 

[ o r + Qj otherwise 
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the set {i,j,k} takes on the value {1,2,3}, {2,3,1} or {3,1,2}, and the same 
is true for the set {/, m,n}, and 
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Abstract 

Three different versions of the Fast Multipole Method (FMM) are employed to 
reduce the storage and computational requirement of the boundary integral in the finite 
element-boundary integral method. By virtue of its 0(A r15 ) or 0(A r1 - 33 ) operation 
count, the application of the FMM, results in substantial speed-up of the boundary 
integral portion of the code, independent of the shape of the B1 contour. The main 
goal of the paper is to provide a comparison of the various FMM approaches on the 
basis of CPU time and accuracy. We present such comparisons and draw conclusions 
on the basis of computing the scattering from large grooves. 
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1 Introduction 


I he finite element -houndary integral ( FE-BI ) method has been quite popular and extensively 
applied to many applications. The method [l] combines the geometrical adaptability and 
material generality of the FEM with the rigorous BJ for truncating the mesh Nevertheless 
although -exact", the FE-BI leads to a partly full and partly sparse system and is thus 
computationally intensive for large boundary integrals. When the boundary is rectangular 
or circular, the FFT can be used to reduce the memory and CPE requirements down to 
A log A [1 ].['-]• However, in general, the boundary integral is not convolutional and in 
that case the CPE requirements will be of 0(.Y 2 ). where \ h denotes the unknowns on the 
boundary. The application of the fast multipole method enables the computation of the 
boundary matrix-vector product with a 0(A' 15 ) or Of.Y 1 33 ) operation count [3]. [4], 

In this paper, we apply three different versions of the Fast Multipole Method (FMM) to 
reduce the storage and computational requirements of the boundary integral when the size 
of the contour is large. By virtue of its low operation count, the application of the FMM 
results in substantial speed-up of the boundary integral portion of the code independent of 
the boundary shape. 


2 Problem Definition and Formulation 


As an application of the proposed technique, we consider the scattering by a cavity- backed 
aperture shown in Figure 1. The FE-BI formulation for this problem was already outlined 
in [2] and results in the solution of the system 


' [A'i m ] r {*} \ _ / o i 

m pi] J i {<m / ~ \ {*«} / 


(i) 


The vector and typical form of this system is given in Figure 2. {<j>} represents the magnetic 
field at the nodes within the groove and on the boundary. Also, {0, } represents the magnetic 
currents on the boundary. By virtue of the finite element method, the matrices [I<] and {B x } 
are sparse and thus the corresponding matrix-vector products are implemented using O(N) 
operations. However [Pi] is a full sub-matrix and thus 0(N *) operations are needed to 
perform the product [Pi]{^} with N b denoting the number of nodes on the cavity aperture. 
Consequently, in an iterative solution, this matrix-vector product becomes the computational 
bottleneck. To reduce the operation count we will herewith employ the FMM procedure for 
implementing the products [Pi]{t/»i}. 

Next we examine three versions of the FMM to accelerate the boundary integral matrix- 
vector product computation. Specifically, the exact FMM [5], [6], a windowed FMM [7] and 
an approximate FMM [8] are examined. 


3 FMM Techniques 

As stated above, the FMM is a fast method for calculating the matrix-vector product. The 
computation of the matrix-vector product is illustrated with the boundary integral for TE 
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Figure 1: Geometry of the groove recessed in a ground plane 

incidence (H-pol). In this case we deal with the discretization of the integral equation 

HI + H' Z = H[ em (2) 

with H s z the scattered magnetic field, H l z the incident magnetic field and H EEA1 the total 
field which is employed to enforce tangential continuity and 

H? = - ^ J rt (*o|? - 7|) if (3) 

M, = E, is the magnetic current and Tj indicates the extent of the boundary, p and p' are 
the vectors describing the observation and source points. 

3.1 Exact FMM 

As illustrated in Figure 3, the Nb boundary unknowns are subdivided into groups with each 
group assigned Mb unknowns. Thus a total of Lb « ^ groups are constructed. Next, by 
invoking the addition theorem for cylindrical wave functions, the Hankel function is expanded 
as 

_ Q/2 _ 

H!> 2) (ko\W + p-p'\)* E U^\p' -p\)H^\k 0 p w )e^'-^'^ (4) 

»=—(?/ 2 

where 4>w and <f> p > p are the angles ~pw and p' — p make with the x-axis respectively and is valid 
for pu> > | p' — ~p\. It should be noted that the source and observation point vectors, p' and p 
have their origin at the center of the source and test groups respectively. The semi-empirical 
formula 

Q/2 = k 0 D + 5ln(k 0 D + tt) (5) 

where D is the diameter of the circle enclosing the groups is used to truncate (4) and in 
general, Q/2 = Mb , assures convergence. The Fourier integral form of the Bessel function 

Jn{ko\7-p\) = ^~ [ ( 6 ) 

Z7T J2ir 
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column 


F i g ure 2: Typical form of the FE-BI system matrix arising from the scattering/radiation 
problem of a groove in a ground plane. 

in conjunction with (4) is used to write (3) as 

H z {p)— J 2 Vi’( < t > )Tii'(<f>)e~ :iko ' p d<t> (7) 

where the far-field pattern of the source group is given by 

= M z {7y r °' ? dl' ( 8 ) 

and the translation operator is given by 

e/2 

T,v{4>) = n^HkoPu (9) 
n=-e/2 

with the propagation vector of the plane wave 

= ko(x cos <t> -f y sin 4>) pQ) 

The integral over 4> is discretized into Q plane wave directions, thus yielding an expression 
for the fields radiated by the source group V to the receiving group / 

H S Ap) = -¥r**Y.T l r(4,,)V e (4,,) e -X* m, 
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where Ao — ~~/Q indicates the angular s|>anng between propagat ion vectors of plane waxes 
emanating from a group and thus ©, = </Ao q = 1 . . . Q and I” = A-„i i cos o . - ( }sin s. 
As mentioned earlier, the number of plane wave directions is set equal to twice the numbei 
of elements in the group ( Q = 2.U,.). thus satisfying the Nyquist sampling theorem with 
respect to the integration over o. 
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Figure 3: Computation of the boundary integral matrix vector product using exact FMM 

In the exact FMM, the matrix vector-product is computed in the following sequence 

1. The pattern of the source group is computed. Mathematically, this corresponds to 

Vi>{<f> q ) = [ M z (p ! )e }kqp ' dl' (12) 

Evaluating this pattern for a single source group and in a single direction requires M b 
operations, corresponding to the number of elements in the group (the integration over 
the line segment is performed as a summation). Consequently for Lb groups and Q 
directions for each group, the operation count is QMbL b . 

2. The translation operation is employed to evaluate the pattern of a source group at a 
test group. Mathematically, 

Mfa) = V v {<f> q )T w {(t> q ) (13) 

For Q directions and Lb source and test groups, this process involves an operation 
count of QLl. 

3. At the receiving group, the fields are redistributed. Mathematically, 

H z (?) = (14) 

q=l 

Thus, computing the redistributed field at a single point requires Q operations, and 
for L b groups each containing M b unknowns the operation count is QL b M b . 
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I he total operation count in the above sequence is (\Q M. I.. - ( ' .Q L; . Hy dioosinc. Q - M. 

for convergence. the operation count is given l.v (\.\L.X, ~C;~. and on sett mo \/. ^ v ~ 

the final operation count is .Y> \ Improvements to the opera, ton count can he achieved bv 
nesting groups leading to the multi-level FMM. 

Beyond the math, the above breakdown of operations are based on the manacer- worker 
model. Basically, we can view each group as managed by the center element with , he workers 
comprising the elements of the group. Communication/interaction among the groups takes 
place through the managers which in turn interact with the group elements. However, this 
type of model is based on certain simplification/decompositions of the original boundary 
integral. Clearly, they reduce the direct interdependence of each group member with other 
elements belonging to different groups. This is the essence of the CPC speed up afforded 

n FMM. However there are inherent approximations which must be understood in order to 
assess the accuracy of each FMM algorithm. 


3.2 Windowed FMM 

In the exact FMM, the translation operation between groups was considered isotropically 
But it is suggestive that the groups would interact very strongly along the line joining them 
and less strongly in other directions. It was shown that by employing a high frequency 
analysis [7], that the translation operator could be contemplated as being composed of a 
geometrical optics term, along the line joining the source and test group, and two uniform 
diffraction terms associated with the shadow boundaries of the GO term. The translation 
operator for different group separation distances along the groove of width 50A is shown in 
igure 4. The number of unknowns on the boundary are 750, resulting in 27 groups. It is 
seen that the lit region of the translation operator decreases with increasing group sepa- 
ration distance, eventually displaying the predictable sine function behavior for large group 
separation distances. The high frequency analysis enables the identification of a lit region 
even for groups which are not very widely separated. Thus the plane wave interactions can 
be filtered by defining a filter function as 


Ww(<b n ) = i 1 (I <t>q ~ (f>u>\ < (3 S ) 


with 


A = sin " (U£) 

and o is the taper factor. The discretized plane wave expansion i 


is now 


Hf = - 


koYo 

4tt 


A<f> £ Wui^Tiri^V 

9=1 


(15) 

(16) 

(17) 


The operation count now associated with (13) is now reduced to C 3 LI ~ with the cor- 
responding total operation count given by C^MbNb + C 4 §j. Grouping the unknowns with 
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A / per group, results in an operation count of Qi A ^ "i. The computation of the bound 
ary integral matrix vector product In employing the windowed I MM is depicted pictorial h 
in Figure 5. It is seen that the filter has the effect of eliminating plane wave interaction: 
at directions away from the line joining the interacting groups. The spectrum of interne 
tions around the line joining the centers of the two groups which are retained, reduces will 
increasing group separation distance. 



0=^4; Pff'— 

0=54; i^= 

Q=54;P^.,= 


<t> - ^(degrees) 


Figure 4: The Translation operator for different groups on the boundary of a 50A wide 
groove; 750 BI unknowns; 27 groups 


3.3 Fast Far Field Algorithm 

The third algorithm employed for hybridization was the fast far field algorithm (FAFFA) 
which is an approximate version of the FMM. Unlike the exact FMM, where the kernel is 
approximated with the addition theorem, in this algorithm the large argument approximation 


H { 0 2) (k o 


P~ P 


e 


-jkoPtU'Plm 


l~2\ e - i k oPiu . ^ _ 

-Plm . I J — J^O PltfPnl* 


7T \/7copiu 


(18) 


is used, where pm is the distance between the center of the test group l and the center of 
the source group p ni , is the distance between the nth source element and its group center 
and pi m is the distance between the mth test element and its group center, as depicted in 
Figure 6. It should be noted that the use of the large argument expansion, as an additional 
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Figure 5: Computation of the boundary integral matrix vector product using windowed 
FMM 

approximation, necessitates that the FMM procedure can now be used for groups which are 
only very well separated. The decoupling of the test-source element interactions in the kernel 
as in (18) enables the computation of the matrix- vector product for far-field groups with a 
reduced operation count as detailed in the following sequence. 

1 . For each test group, the aggregation of source elements in a single source group involves 
M b operations, corresponding to the number of elements in the source group. The 
aggregation operation corresponds to 


V vl = £ Mje~* 




2. Since the above aggregation operation needs to be done for all source groups the 
operation count becomes 0{^M b ) ~ 0{N b ), where ^ represents the total number of 
groups. Also this operation, being dependent only on the test group rather than the 
test element, needs to be repeated for ^ test groups leading to a total operation count 

of ^(a^) for aggregation. 

3. The next step is a translation operation corresponding to 


where 


Aft = Ti'iVf 


2j 


T ltl = s^ 


7T VW/'Z 


Since this needs to be done only at the group level, it involves operations for all 

possible test and source group combinations and is the least computationally intensive 
step. 
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Figure 6: Computation of the boundary integral matrix vector product using the FAFFA 


4. The final step in the sequence would be the process of disaggregation corresponding to 
the operation 


o 


Nb/Mb 


H 1{P) = 7T~ 53 Ane j *° p,Vp '" 


/'=1 


( 22 ) 


Conceptually, this process is the converse of aggregation. Since this operation involves 
only the source group instead of the source element it needs to be done for each source 
group thus implying an O(^) operation to generate a single row of the matrix-vector 
product. To generate A/(, rows corresponding to a test group the operation count would 
be O(Nb). With ^ test groups, the operation count would be of 0{jf). 


5. The near field operation count being of 0{N b M b ) and the far field being gives 

a total operation count of 6 


Op. count ~ CiN b M b + C 2 ~ 

M b 


(23) 


Typically, we can set the elements in each group, M b — y/ N b and as a result the total 
operation count is 0 


4 Results 

A computer code based on the above formulation was implemented and executed on a HP 
9000/750 workstation with a peak flop rate of 23.7 MFLOPS. The geometry considered was 
the rectangular groove shown in Figure 1. The performance of the hybrid algorithms with 
respect to accuracy and speed were compared. The benchmark for accuracy was the RMS 
error [9]. Table 1 compares the execution time and error of the standard FE-BI to the 
FE-Exact FMM, FE-FAFFA and the FE-Windowed FMM for grooves of widths 25A, 35A 
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and it) A. 1 lie depth of the groove was 0.35A with a material fillmc of ,. = 1 and ,i =1 
and was illuminated at normal incidence. The data reveals that the FK-Kxact FMM 

rr' 1 ?! rr? 1 ' 7 SavingS in execution timr V1,h compromise in accuracv. While the 

L-I AFFA is the fastest of the three algorithms, the RMS error was substantially hi-dior 
(>1 dB). If the maximum tolerable RMS error is set at ] dB [O’. the FF-Wmdowod FMM 
is the most attractive option as it meets the error criterion while being onlv sliehtlv slower 
than the FE-FAFFA. An observation of interest was the comp arisen of the residual error as 


Groove 

Width 

Total 

Unknowns 

BI 

Unknowns 

^ ^ [ T tint for BI ( Minute s, st couds) I 

FE-BI (C’G) 

FE-FAFFA 

FE-Exact FMM 

FE-WFM.M 

25 A 

2G31 

375 

(8.48) 

(3.26) 

(5.25) 

(4.13) 

35A 

3681 

525 

(16.34) 

( 5.55 ) 

(10.31) 

(7.22) 

50 A 

5256 

750 

(45.1) 

(14.31) 

(26.18) 1 

(16.10) 



RMS error (dB) 1 

Groove Width 

FE-FAFFA 

FE-Exact FMM 

FE-WFMM 

25 A 

1.12 

0.0752 

0.6218 1 

35 A 

1.2 1 

0.1058 

0.721 

50A 

1.36 

0.1123 

0.843 


Table 1: CPU Times and RMS error of the hybrid algorithms 

a function of the^number of iterations in the CG solver for the hybrid algorithms. This is 
shown in Figure 7. It is seen that the convergence curves for the FE-BI and the FE-Exact 
FMM overlap each other to graphical accuracy while the FE-Windowed FMM shows a very 

small deviation. Thus, the hybridization of the FMM does not have a deleterious effect on 
the conditioning of the FE-BI system. 

The execution time of the fastest of the three hybrid techniques, the FE-FAFFA is com- 

Pa m d u7 it o h t ! ie FE ' BI em P lo y in § the s P e dal CG-FFT solver, suited for only planar apertures 
in Table 2. It is seen that the CG-FFT solver is substantially faster but is applicable only 
to convolutional boundary integrals. 

The performance of the hybrid algorithms at a more stressing angle of incidence is de- 
picted in Figure 8. The RMS error follows the same trend as for normal incidence illumi- 
nation. The width of the groove illuminated was 10A and this example reveals that the 
technique is scalable for smaller problems. The near group radius in this example was 1A- 
implying that the matrix vector products for groups separated by a distance less than a 
wavelength was computed using the exact method of moments procedure. It has to be noted 
that the application of the hybrid techniques for the 10A groove illustrates that the near- 
group radius can be reduced as the problem size becomes smaller down to a minimum in the 
vicinity of 0.3A. 
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Figure 7: Convergence curves for the hybrid algorithms for the groove of width 25 A 

4.1 Summary 

The hybridization of the finite element - boundary integral and fast multipole methods was 
examined from a computational performance point of view. Three different versions of the 
fast multipole method was used for executing the matrix-vector product associated with 
the boundary integral. It was shown that a considerable reduction in CPU time could be 
achieved and further reductions are likely as the boundary integral increases in size. The 
FE-Windowed FMM provides the best compromise in terms of speed and accuracy. The 
FE-BI with the CG-FFT solver is faster than any of the FEM-FMM versions and can be 
used if the terminating boundary is amenable to its use. 
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i ( pi } 

' inn for HI . Muitih 
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i I L-Bl iCG 

i FE-FAFFA IT 

i HI i CGFFT i 

20 A 
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O.-lSi 

1 3.26 i 

-1.41 - 

30 A 

36S1 
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IK..31) 
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(c) 


Figure 8: Scalability of the hybrid techniques to smaller problems (a) Problem geometry (b) 
Bistatic patterns (c) Error table 
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Abstract 

A metal-backed layer of absorbing material offers a number of advantages for truncat- 
ing the computational domain in a finite element simulation. In this paper we examine 
isotropic and anisotropic absorber layers for the purpose of truncating finite element 
meshes. Optimal design curves are presented for these absorbers that can be used to 
select the various parameters (thickness, propagation and sampling rate) so the reflectiv- 
ity is minimized. Applications to radiation and scattering problems are also considered, 
and these illustrate the accuracy and versatility of the absorber layers for large scale 
electromagnetic simulations. 


1 Introduction 

In the numerical solution of electromagnetic scattering and radiation problems it is necessary' 
to truncate the computational domain in a manner which ensures that the waves are outgoing. 
This is true also in the analysis for many microwave circuits, and the need to terminate the 
mesh is common to finite element (FEM) and finite difference-time domain (FDTD) methods. 
One way to do so is to enforce an absorbing boundary condition (ABC) at a surface as close as 
possible to the scatterer or radiator, and a review of available ABCs has been given by r Senior 
and Volakis [1], Unfortunately, ABCs are limited in their ability to conform to the surface of 
the scatterer. They may also require an a priori knowledge of the wave’s properties and, in an 
FEM solution, they generally reduce the rate of convergence. Another way of terminating the 
mesh is to use a metal-backed layer of isotropic or anisotropic absorbing material [2, 3, 4, 5], and 
such layers are often referred to as artificial absorbers since their material parameters may be 
physically unrealizable. 

The implementation of artificial absorbers for finite element mesh truncation is illustrated 
in Fig. 2 and as expected the layer’s material composition plays a major role in the perfor- 
mance of the artificial absorber. However, the chosen numerical discretization of the absorber 
has an equally important role and cannot therefore be ignored in the design of the absorber for 
numerical simulations. In this paper, we examine the performance and design of both isotropic 
and anisotropic homogeneous absorbers for truncating finite element meshes. Their applica- 
tion to the finite element solution of three dimensional radiation and scattering problems is 
also considered and results are shown which demonstrate the utility of these mesh truncation 
schemes. 
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2 Analytical Study 

Consider the metal-backed planar layer shown in Fig. 2(a). The surface ,r = 0 is the interface 
between free space (in r < 0) and a lossy material (in x > 0) backed by a perfect electric 
conductor at x = i. For an incident plane wave 

£. or H; = ( — ( J COS C+ 1/ Sill O I J I J 

the reflection coefficient is Re(o) or /?//(©). and the objective is to minimize these. 

If the layer is composed of a homogeneous isotropic material whose relative permit tivitv < r 
and relative permeability p r are such that c r = p r = b = o - j3 (say), then 


Re{<P) = 


Rh(4>) = 


\/l — b 2 sin 2 0 — j cos 0 tan ( k 0 b1\J 1 — 6“ 2 sin 2 0 ) 
yl — h -2 sin 2 0 + j cos <p tan ( kybf yj 1 — 6~ 2 sin 2 0) 
yl — f> -2 sin 2 0 - j cos0cot(/eo&fy^l — 6 — 2 sin 2 </>) 
yl - f> -2 sin 2 0 + j cos 0cot(fc o &<yi — 6“ 2 sin 2 0) 


( 2 ) 


( 3 ) 


These differ because the presence of the metal backing has destroyed duality, and at grazing 
incidence (0 = tt/2), R e = R h = -1. If sin0 > 1, i.e. 0 = 7r/2 + with £ > 0 so that 
sin<£ = cosh 8 and cos 0 = —j sin 5, \R\ differs from unity by only a small amount for both 
polarizations. The behavior of \Re,h(^)\ as a function of sin0 is illustrated in Fig. 3. At 
normal incidence (0 = 0), (2) and (3) give 


Re,h( 0) = =Fe 2 tk 0 t(a-jP) ( 4 ) 

whose magnitudes are independent of a and can be made as small as desired by choosing 
k 0 tfi sufficiently large. Since large values of 0 produce rapid field changes in the dielectric, a 
disadvantage is that high (and often very high) sampling rates are necessary to simulate them 
numerically. 

As k 0 t — > 00 , Re,h{4>) — t 0 only for normal incidence, but Sacks et al [4] have shown that 
a particular uniaxial anisotropic material has this property for all 0 < tt/2. The result is an 
example of a perfectly matched layer (PML), and if 


0 


H r = bl — (b — -)xx 
b 


( 5 ) 


where I is the identity tensor and b = a — j/3 , then 

R E jt{4>) = 

which reduces to (4) in the particular case of normal incidence. If k o 0 » 1 the reflection 
coefficients decay exponentially for all 0 < 7 t/ 2, and since (6) is also valid for sin 0 > 1, the 
choice a > 0 ensures an exponential decay for these angles as well. The behavior of \Re,h{4>)\ 
is illustrated in Fig. 3 for the same values of k 0 t, a and 0 used for the isotropic layer. Clearly, 
a major advantage of the PML is that its reflection coefficient remains low for a wide range of 
angles of incidence. 
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3 Absorber Design 

71k- theoretical behavior of the reflection coefficient |/?j for the metal-backed layers is relatively 
simple. In the case of the isotropic material, an increase in 3 and/or / decreases ( /* > ( 0 ) | . The 
uniaxial material has this behavior for all real angles of incidence, and while o plays little or 
no role, large values of a do produce higher absorption for complex angles. It follows that 
for a uniaxial layer of given thickness /. a and 3 can be chosen sufficiently large to produce 
high absorption over a broad angular spectrum, with angles near o = rr/2 providing the only 
exception. 

Imfoitunateh. the analytical results do not immediately translate into numerical perfor- 
mance. Because of the discretization inherent in an FEM implementation, the fields inside the 
layer are reproduced only approximately, and this is particularly true for a rapidly decaying 
field. To design a good absorber it is necessary to understand the impact of the sampling rate 
on the choice of a, 3 and f, and our objective is to find the minimum number of sampling ele- 
ments ( or discrete layers) to achieve a specified \R\. It is anticipated that the errors introduced 
by the discretization will have a number of consequences. In particular, for a given number N 
of discrete layers and given f, increasing 0 will ultimately lead to an increase in \R\ because 
of the inability to model the increasing attenuation, and an increase in a will likely produce a 
similar effect. To obtain some insight into the roles played by N, a, (3 and t , a simple FEM 
model for computing the reflection coefficient of a metal-backed absorber layer (see Fig. 1 ) was 
considered. 

We examine first a homogeneous isotropic layer at normal incidence for which the theoretical 
reflection coefficients are given in (4). In spite of the fact that the magnitudes are the same 
for both polarizations, a polarization dependence shows up in the FEM implementation. This 
is illustrated in Fig. 4, and we note that as N increases, the FEM values of |i?(0)| converge 
to the common theoretical value for both polarizations. The nulls associated with the H- 
polarization curves are characteristic of the behavior of the modeled absorber layer, and are 
due to the interference of the reflected fields from the dielectric interface, the metal-backing 
and the individual layers used to model the absorber numerically. 

Given the many parameters (/?, a, t, N), it is essential to consider some optimization of 
the proposed metal-backed absorber as a function of these parameters. For this purpose, a 
study was carried out using the FEM code mentioned above. Initially, an investigation was 
performed to examine the effect of each parameter on the absorber’s performance. As can be 
expected, o plays little role in the performance of the absorber for relatively small values of 
Also, larger values of N for the same thickness t lead to lower reflection coefficients. 
An optimum value exists but this depends on /?, the absorber’s decay parameter. However, the 
most important observation from this preliminary study concerns the scalability of the product 
fit. It was found that for small values of a and for thicknesses up to at least 0.5 wavelenths, the 
reflection coefficient is a function of the product (3tj\ 0 alone. As a result, one can construct 
reflection coefficient curves as a function of (it which are optimum for a given N . Fig. 5 gives 
such design curves by plotting |i?(0)| and N versus (3t/ \ 0 on the same figure, and these can be 
used to determine the optimum N for a given (it/X 0 . To see how to use the figure, suppose that 
the desired reflection coefficient at normal incidence is -50 dB. In Fig. 5 we observe that the 
|fl(0)| curve intersects the -50 dB line at (3ij\ 0 ~ 0.58, and referring now to the N curve, the 
number of elements required is A r = 10. The value of (3 can then be found by specifying either 
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thr element size or the layer thickness. 1 hus. for elements 0.025A,. thick, we have t = ().25A t , 
and 3 = 2.32. By increasing A we can improve the performance up to the limit provided bv 

the theoretical value of |/?(0)| which has been included in Fig. '>. A good approximation to the 

long-dashed curves in Fig. 5 obtained by linear regression is 

3t 

— = — 0.01 06| /?j + 0.0433 (7) 

A' = 0.147 exp [7.353J//A 0 ] (S) 

where |/?| is measured in dB and A' is the smallest integer equal to or exceeding the right hand 
side of (8). 


4 Application to Antennas and Scattering 

The above study dealt with planar absorber layers but the layers can also be curved for ter- 
minating finite element meshes in a conformal manner. This is illustrated in Fig. 1 where the 
artificial absorber layer is used to terminate the computational domain surrounding a patch 
antenna. Conformal mesh terminations are quite attractive in FEM modeling because they lead 
to a substantial reduction of the computational volume and are also compatible with structured 
as well as unstructured meshes. In contrast to absorbing boundary conditions, they do not re- 
quire a knowledge of the closure’s principal radii of curvature and uniaxial artificial absorbers 
offer the possibility of reflectivity control as low as -60dB. 

To illustrate the applicability of the artificial absorbers in FEM modeling, two examples are 
considered, one dealing with antenna analysis and the other with scattering by a non-canonical 
structure. For simplicity both cases employed a simple version of a homogeneous artificial 
absorber consisting of three layers 0.05A o thick as shown in Fig. 1. The attenuation constant 
for each layer was 0 = 2.6, and thus g = 0.405. From Fig. 5 this absorber provides a normal 
incidence reflection coefficient of -30dB, and from the same figure we also read off that this 
reduction can be realized with 3 layers or more. Thus, the design of the proposed homogeneous 
absorber is consistent with the curves given in Fig. 5. 

The absorber termination shown in Fig. 1 has been used to model a variety of patch 
(circular and rectangular) antennas on doubly curved platforms, including circular, spherical, 
conical and ogival surfaces. Of particular interest is the computation of the resonant frequency 
which is a rather sensitive quantity and its accurate computation via the proposed FEM model 
provides a good test of the absorber’s performance. Fig. 6 shows the results for a rectangular 
patch antenna mounted on the conical surface illustrated in the figure. The patch resides on 
a substrate having e r = 2.32 and a thickness h = 0.114cm. Its dimensions are given in Fig. 6 
and on the basis of the approximate cavity model it resonates at ZflGHz. From the computed 
input impedance plot, it is seen that the resonance frequency predicted by the FEM code is 
3.115 GHz, which is within 3 percent of the cavity model. The FEM computations were carried 
out using prismatic edge elements and the surface grid is also shown in Fig. 6. A total of 2358 
prisms were used for this analysis resulting in 3790 degrees of freedom. 

As a scattering example we consider the radar cross section of the NASA metallic almond 
[6] shown in Fig. 7. This body is 9.936 inches long, and precise formulae for describing its 
surface are given by Woo et. al . [6] . Measured data are also available at several frequencies and 
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can be used to benchmark the accuracy of the simulation. Our finite element code F EM ATS 
[i] was used to model the almond illuminated with a plane wave at a frequency 1.19(lHz. 
This code employs edge-based tetrahedral elements and the mesh was again terminated using 
the aforementioned 3-layer homogeneous artificial absorber. For ease of mesh generation, a 
structured prismatic mesh was first generated conformal to the almond s surface and consisted 
of nine 0.05A layers with the outer layers occupied by the artificial absorber. The structured 
prismatic mesh was then turned into a tertrahedral mesh resulting in a total of 40.878 edges. 
Fig. 7 displays the computed radar cross section(RCS) computations with the measured ones 
for both polarizations of incidence. The patterns are taken in the plane most parallel to the flat 
side of the almond(non-s\ mmetnc plane), with zero degrees corresponding to incidence tip-on. 
As seen, the calculations are in good agreement with the measured data except at near 90 
degrees for HH polarization. However, other reference calculations based on a moment method 
code are in agreement with the FEMATS data, suggesting that the discrepancy may be due to 
minor alignment errors in the measurement. 

5 Conclusion 

In this paper we investigated homogeneous isotropic and uniaxial artificial absorbers for finite 
element mesh truncation. By properly selecting the material properties and sampling rate, it 
was demonstrated that almost any desired level of absorption can be attained, and tvpicallv 
very few samples (less than five) are needed to achieve a reflection coefficient of -40 dB over 
a wide range of incidence angles. Design curves were presented which can be used to select 
the various parameters (loss, thickness and sampling rate) on the basis of a desired reflection 
coefficient. As expected, a lower material loss requires a thicker absorber to produce the same 
reflection coefficient but, on the other hand, a higher attenuation rate requires more samples to 
attain a lower reflection coefficient in a numerical implementation. Most likely, inhomogeneous 
(tapered) artificial absorbers can lead to lower reflection coefficients, but these have not yet 
been investigated to any great extent. 

In contrast to absorbing boundary conditions, a particular advantage of the proposed ab- 
sorbers is that they can be used to terminate a finite element mesh conformal to the target 
or radiator surface without needing a priori information about the wave’s propagation char- 
acteristics. To test the performance and applicability of the proposed absorber for truncating 
finite element meshes in a three dimensional setting, two examples were considered. Namely, 
computation were carried for the input impedance of a patch antenna on a conical surface and 
the radar cross section of a non-canonical slender body. In both cases the computed values 
were in good agreement with reference data by using the proposed artificial absorber for mesh 
truncation. 
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Figure 2: Geometry of the metal backed absorber layer and its finite element discretization. 



Figure 3: Analytical results for homogeneous isotropic and anisotropic absorbing layers with 
t = 0.25A o and 6=1 — j 2: (• • •) isotropic E pol, (- • -) isotropic H pol. and ( ) anisotropic. 
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Figure 4: Numerical results for a homogeneous isotropic and anisotropic absorbing layer with 

t = 0.15A o and b = -j 2.5. The top four curves are for Epol., the bottom four for H pol • ( ) 

exact, (- - -) N=3, ( ) N=6, and ( ) N=12. 
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Figure 5: Absorber design curves. The straight lines give |fl| in dB, and the curved ones give 
N '■ ( ) exact |/?|; ( ) numerical |/?| 
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Figure 6: Microstrip sectoral patch on cone: (a) configuration, (b) patch input resistance as a 
function of frequency for different cone angles and (c)display of the mesh vertical fields below 

the natch and in the substrate 
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RCS - NASA METALLIC ALMOND 



(b) 


Figure 7: Geometry and RCS results for the NASA almond[6!. (a) surface and volume prismatic 

conformal mesh of the 9.936 almond and (b) RCS pattern calculations and measurements at 
1.19GHz. 
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ABSTRACT 

A metal-backed layer of absorbing material offers a number of ad- 
vantages for truncating the computational domain in a finite element 
simulation. In this paper we present design curves for the optimal se- 
lection of the parameters of the layer to achieve a specified reflection 
coefficient. The curves are based on one-dimensional finite element sim- 
ulations of the absorbers, and the optimization is therefore a function 
of the sampling rate. Three types of material are considered, including 
the recently introduced perfectly matched uniaxial material, either ho- 
mogeneous or with a quadratic material profile. Two three-dimensional 
applications are also presented and used to examine the validity of the 
design curves. 

1 INTRODUCTION 

In the numerical solution of electromagnetic scattering and radiation 
problems it is necessary to truncate the computational domain in a 
manner which ensures that the waves are outgoing. This is true also 
in the analysis for many microwave circuits, and the need to terminate 
the mesh is common to finite element (FEM) and finite difference- 
time domain (FDTD) methods. One way to do so is to enforce an 
absorbing boundary condition (ABC) at a surface as close as possible to 
the scatterer or radiator, and a review of available ABCs has been given 
by Senior and Volakis [1]. Another way is to use a metal-backed layer 
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of isotropic absorbing material J.-i . but both sc hemes have limitations, 
lor example, an AB( recpiires a prion knowledge of the propagation 
constant which, in a microwave problem, may not be the same in all 
section of the computational domain. Also, when used to terminate 
an open domain. ABCs reduce the convergence rate and may be hard 
to implement on a surface conformal to the scatterer or radiator. An 
isotropic dielectric layer alleviates some of these difficulties, but its 
accuracy and aspect coverage are limited. 

Recentlc a new anisotropic absorber has been proposed for termi- 
nating the domain. By introducing an additional degree of freedom. 
Sacks et al. [4] have shown that a uniaxial material can be designed to 
have zero reflection coefficient at its interface for all angles of incidence. 
If the material is also loss}', a thin metal-backed laver can be used to 
terminate an FEM mesh, and though the material is no longer realiz- 
able physically, the associated fields are still Maxwellian. This is often 
referred to as a perfectly matched layer (PML). and its development 
was motivated by the non-Maxwellian layer introduced by Berenger [5] 
(see also [6]) for FDTD problems. By choosing the parameters ap 
propriately, it is possible to achieve any desired level of absorption for 
almost all angles of incidence using only a thin layer, but its numerical 
simulation is a more challenging task. Because of the rapid exponential 
decay of the fields within the layer, there are large variations in a small 
distance, and it is difficult to reproduce these in a numerical simula- 
tion. Thus, for a discretized PML, the numerical sampling as well as 
the material properties affect the reflection coefficient that is achieved. 

In this paper we consider the design and performance of three types 
of metal-backed planar layers for terminating FEM meshes: homoge- 
neous isotropic, homogeneous anisotropic (uniaxial), and inhomoge- 
neous (tapered) uniaxial materials. Using one-dimensional finite ele- 
ment simulations, their numerical performance is examined and com- 
pared with their theoretical capability. Not surprisingly, the sampling 
rate has a major effect on the reflection coefficient. Based on a detailed 
numerical study, we identify scalable parameters in the numerical model 
and use these to generate design curves and formulas for choosing the 
sampling rate and material properties to achieve a specified reflection 
coefficient. As expected, a tapered uniaxial material proves superior 
to the homogeneous one. The applicability of these results to three- 
dimensional problems is then illustrated for a simple microwave circuit 
and a rectangular waveguide. 
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(a) 


Figure 1: Geometry of the metal-backed absorber layer (a) and its FEM 
implementation (b). 

2 ANALYTICAL STUDY 

Consider the metal-backed planar layer shown in Fig. 1(a). The surface 
x = 0 is the interface between free space (in x < 0) and a lossy material 
(in x > 0) backed by a PEC at x = t. For an incident plane wave 

E. or H~ = e -i k o(x COS <t>+y sin 4 i) 

the reflection coefficient is Re(4>) or Rh{4>), and the objective is to 
minimize these. If the layer is composed of a homogeneous isotropic 
material whose relative permittivity t r and relative permeability fi T are 
such that e r = = b = a — j/3 (say) where a and /? are real, then 



These differ because the presence of the PEC backing has destroyed 
duality, and at grazing incidence (4> = 7t/2), Re = Rh = — 1. If 
sin^ > 1, i.e. (f> = ~/2 + jS with 8 > 0 so that sin d> = cosh 8 and 
costfi = —jsm8, |7?j differs from unity by only a small amount for 
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sin 0 

Figure 2: Analytical results for homogeneous isotropic and anisotropic 
absorbing layers with t = 0.25A o and b = 1 - j2: ( ■ - •) isotropic E pol., 
(- • -) isotropic H pol. and ( ) anisotropic. 

both polarizations. The behavior of |/?e,//( 0)| as a function of sin 0 is 
illustrated in Fig. 2. At normal incidence (0 = 0), (2) and (3) give 

Re.h{ 0) = ^ e ~ 23kot{a - J0) (4) 

whose magnitudes are independent of a and can be made as small as 
desired by choosing k 0 t/3 sufficiently large. 

As k 0 t -> oo, Re,h{4>) — > 0 only for normal incidence, but Sacks et. 
al. [4] have shown that a particular uniaxial anisotropic material has 
this property for all 0 < tt/ 2. The result is an example of a perfectly 
matched layer (PML), and if 

l r = f T =bl-{b- l -)xi (5) 

where / is the identity tensor, then 

Re,h{ 4>) = ^ e ~ 2jkot(a ~ }0)cos,p (6) 

which reduces to (4) in the particular case of normal incidence. If 
ko(3 1 the reflection coefficients decay exponentially for all 0 < 7r/2, 
and since (6) is also valid for sin0 > 1, the choice o? > 0 ensures an 
exponential decay for these angles as well. The behavior of \R E} H{(f>)\ 
is illustrated in Figure 2 for the same values of k 0 t , a and (3 used for 
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the isotropic layer. Clearly, a major advantage* of tin* 1W1L i> that its 
i e fleet ion coefficient remains low for a wide ranue ol ancli's of incidence. 

Although the outer surface of t lie layer is reflect ioniess for all o. 
the abrupt change in the* material properties at .r = 0 may produce a 
contribution in an FEM solution. We can eliminate t Ik* discont iiiuit v i>v 
tapering tlie properties as a function of x to produce an inhomogeneous 
anisotropic layer. As shown by Legault and Senior [ 7 ], if h = 
a wave propagating into the material has the form 

c -.M'o {^(r) cos o+y sin 0} ^ j 

and when 


l ( *) — 1 + (o — ji 3 — 1 ) j . (S) 

which tends to unity as x — > 0+, the reflection coefficients of the laver 
are identical to those given in (6). With this expression for 7(0*), the 
attenuation is less where the field is larger, i.e. close to the interface, 
and increases as the field is absorbed. A simplified version of (8) is 
employed in Section 3.4. 


3 NUMERICAL STUDY 

For all three types of layer the theoretical behavior of |/?| is relatively 
simple. In the case of the isotropic material, an increase in j3 and/or 
t decreases |/?(0)|. The uniaxial material has this behavior for all real 
angles of incidence, and while a plays little or no role, large values of 
a do produce higher absorption for complex angles. It follows that for 
a uniaxial layer of given thickness t, a and (3 can be chosen sufficiently 
large to produce high absorption over a broad angular spectrum, with 
angles near </> = 7 t/ 2 providing the only exception. 

Unfortunately, the analytical results do not immediately translate 
into numerical performance. Because of the discretization inherent in 
an FEM implementation, the fields inside the layer are reproduced only 
approximately, and this is particularly true for a rapidly decaying field. 
To design a good absorber it is necessary to understand the impact of 
the sampling rate on the choice of a, (3 and f, and our objective is to 
find the minimum number of sampling elements (or discrete layers) to 
achieve a specified |fl|. It is anticipated that the errors introduced by 
the discretization will have a number of consequences. In particular, 
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Figure 3: Numerical results for a homogeneous isotropic layer with 
t = 0.15A o and 6 = —j 2.5. The top four curves are for E pol., the 

bottom four for H pol.: ( ) exact, ( ) N=3. ( ) N=6 and ( — 

— ) N=12. 


for a given number N of discrete layers and given t, increasing / 3 will 
ultimately lead to an increase in |/?| because of the inability to model 
the increasing attenuation, and an increase in o will likelv produce a 
similar effect. To obtain some insight into the roles played by A r , o, (3 
and t, we now consider a simple FEM model of the layers. 


3.1 Numerical Model 

A one-dimensional FEM code was used to examine the numerical per- 
formance of the absorbing layers. The computational domain was lim- 
ited to the discretized layer structure shown in Fig. 1(b), with the 
appropriate boundary conditions applied at the interface x = 0 and the 
PEC backing x — t. 

We consider first a homogeneous isotropic layer at normal incidence 
for which the the theoretical reflection coefficients are given in (4). In 
spite of the fact that the magnitudes are the same for both polariza- 
tions, a polarization dependence shows up in the FEM implementation. 
This is illustrated in Fig. 3, and we note that as N increases, the FEM 
values of |/?(0)| converge to the common theoretical value for both po- 
larizations. 
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3.2 Dependence on o and 3 

Tor a layer of constant thickness the theoretical value of ,/i’M.M is in- 
dependent of a and polarization, but in the numerical implementation 
the behavior is much more complicated. Figure I shows | /,’( (» 1 1 plotted 
\ ersus o and 3 for a layer of thickness / = 0.25A<j made up of 5 ( = .Y) 
elements, where the darker tones indicate lower values. For small A the 
results are in close agreement with theory. As evident from the level 
lines. | /?( 0 ) | is almost independent of o and decreases exponentially 
with 3. leading to a linear decrease on a dB scale. For large 3. how- 
ever. the behavior is quite different, and the most striking feature is the 
series of deep minima whose spacing in o increases with increasing o 
and decreasing 0. These are numerical artifacts which are common to 
both polarizations and may depend on the particular numerical code 
employed. The minima for the two polarization are interlaced, and for 
H polarization the first minimum occurs at o = 0. 3 = 1.6. Their 
locations also depend on t and N. If X is fixed, the spatial sampling 
is inversely proportional to t. Decreasing / results in better sampling, 
pushing the minima to higher values of 0 and producing agreement 
with the theoretical values for larger 0 than before. Increasing t has 
tin opposite effect. On the other hand, if t is fixed, increasing N im- 
proves the accuracy, and shifts the minima to higher 0. Apart from 
the minima, the reflection coefficients for fixed 0 increase slightly with 
increasing a, and it is therefore sufficient to confine attention to the 
lower values of a. 

In Figure 5 the reflection coefficients are plotted as functions of 0 
for the same layer with a = 0 and a = 0.75. The curves correspond 
to vertical cuts through the patterns in figure 4, and we also show the 
theoretical value obtained from (4). We observe that as 0 increases the 
reflection coefficients decrease initially at almost the same rate implied 
by (4), but beyond a certain point they begin to increase. The deep 
minimum at a = 0 and 0 = 1.6 in Figure 4(a) is clearly seen, but 
for design purposes it is logical to focus on the worst case, i.e. the 
polarization for which the reflection coefficient is larger. The upper 
curves in Figure 5 are almost identical and constitute this case. Since 
they correspond to two different values of a, either of them would 
suffice, but for reasons that will become evident later, we choose q = 0. 
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Figure 4: Plot of |/2(0)| in dB for (a) an H polarized and (b) an E 
polarized wave incident on a homogeneous isotropic layer with t = 
0.25A o and N = 5. The solid curves are level lines. 






Figure 5: |/?(0)| with t = 0.25A o and N = 5: ( ) Exact. ( ) 

o = O^E pol . , ( ) q = 0 H pol., (- - - ) q = 0.75 E pol. and (- • •) 

q = 0.75 H pol. 

3.3 Dependence on /?, N and t 

We now seek a connection between the values of 0, A' and t for which 
|/?(0) | is minimized. To this end, we first examine |tf(0)| as a function 
of 0 and N for constant t, and the resulting plot is shown in Figure 6 
for E polarization with t = 0.25A o and ct = 0 as before. For fixed 0 
the reflection coefficient tends to its exact values as N increases. This 
is evident from the level curves and, as expected, the convergence is 
better for the smaller 0. Consider now the behavior of |fi(0)| for fixed 
N . As 0 increases from zero, the reflection coefficient decreases to a 
minimum and then increases. The location of the minima is indicated 
b> the solid line. This is consistent with the behavior shown in Figure 
5 and the upper curve is, in fact, just a vertical cut through Figure 6 
at N = 5. The solid line in Figure 6 therefore gives the value of 0 at 
which |/?(0)| is a minimum as a function of the number of elements. 

If the process is repeated for other layer thicknesses, it is found that 
for minimum |/?(0)| the curve of 0t/X o versus N is virtually the same 
for all thin layers. The observation that 0tj Ao is a scalable parameter 
is an important conclusion of our study, and by choosing a constant 
layer thickness we can produce a universal curve for the optimal choice 
of N and 0 in FEM simulations. Such a curve is shown in Figure 7 
and can be interpreted as giving the value of 0t/\ o for a specified N 
to minimize the reflection coefficient |/2(0)|. For example, if t = 0.2A o 
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0 


Figure 6: |/?| as a function of 0 and N for E pol. with t = 0.25A o and 
Q = 0. The dashed curves are level lines and the solid curve gives the 
location of the minimum |/?|. 

and N = 3, then 0 = 2.13. If a smaller value of 0 is chosen, \R{0)\ 
will be larger (see Figure 5), and this can be attributed to the fact 
that the field reflected from the metal backing has not been attenuated 
sufficiently. If 0 is set to a value larger than 2.13, |/?(0)| will still be 
larger because the chosen N is too small to reproduce the rapid field 
decay within the layer. 

So far we have considered only normal incidence, but for the anisotropic 
layer it is a simple matter to extend the results to all real angles of in- 
cidence <p < 90°. As evident from the exponent in (6), the absorption 
at normal incidence is reproduced at any angle if the layer thickness is 
inversely proportional to cos <f>. This can be achieved by specifying the 
layer thickness t as a fraction 8 of the wavelength A x along the normal 
(or x axis) to the material interface. Since t = 8X 0 /cos(j) = 8X x , t/X x 
is now independent of <f> and the scalable parameter 0tj A 0 (at normal 
incidence) becomes 0tjX x . For the anisotropic layer, all the results 
obtained at normal incidence are made applicable for arbitrary by 
substituting A x for Ao- For example, plotting the optimum 0tjX x as a 
function of N duplicates the curves shown in Figure 7. This notion can 
also be used to account for problems where the outer medium is not 
free space. In such cases, we have A x = A 0 />/c7 (at normal incidence) 
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Figure 7: 3t/X 0 computed from the E pol. case with o = 0: ( ) 

t = 0.15A o . ( ) t = 0.25A o and ( ) t — 0.5A o . 


where c T is the relative permittivity of the outer medium. 

Although the scaling property of (3t/ X x has been established only 
for q = 0. it holds to a reasonable degree for small o / 0, but as a 
increases, the 3tj X x versus N curves become increasingly dependent on 
a. The scalability also extends to the associated values of |/?(<p)|, and 
this enables us to provide a simple design prescription for an absorbing 
layer. 

3.4 Design Curves 

Since the quantities / 3t/X x and |i?(<£)| are the same for layer thicknesses 
up to about 0.5A X at least, design curves can be obtained by plotting 
I R{4 >) I and N versus (3t/X x on the same figure as shown in Figure 8. To 
see how to use the figure, suppose that the desired reflection coefficient 
at normal incidence is -50 dB. At (f> = 0 we now have X x = A 0 . In 
Figure 8 we observe that the |i?(0)| curve intersects the -50 dB line 
at 0t/X o ~ 0.58, and referring now to the N curve, the number of 
elements required is N = 10. The value of (3 can then be found by 
specifying either the element size or the layer thickness. Thus, for 
elements 0.025A o thick, we have t = 0.25A o and (3 = 2.32. By increasing 
N we can improve the performance up to the limit provided by the 
theoretical value of |i?(0)| which has been included in Figure 8. A good 
approximation to the short-dashed curves in Figure 8 obtained by linear 
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Figure 8: Absorber design curves. The straight lines give \R\ in dB. 

and the curved ones give A: ( ) exact |/?|, ( ) homogeneous case 

and ( ) inhomogeneous case. 


regression is 


y = -0.0106|fl| +0.0433 (9) 

N = 0.147 exp [7.353/fr/A*] (10) 

where A x = A 0 / cos 4>, |/2| is measured in dB and N is the smallest inte- 
ger equal to or exceeding the right hand side of (10). These equations 
hold for 4>, 0° < 4> < 90°, for the anisotropic layer, and <f> = 0° for the 
isotropic one. The design criterion provided above applies to specific 
angles of incidence. In the case where a specific absorption level is re- 
quired over a range of angles of incidence the layer must be designed 
for the largest angle occurring. Doing so ensures that the absorption is 
superior at all smaller angles. 

The performance can be improved by making the anisotropic mate- 
rial inhomogeneous, and to illustrate this we consider the case *y{x) = 
— J/3 {x/t } for which the theoretical reflection coefficient is the same as 
before. The scalability is still preserved and the resulting curve is shown 
in Figure 8. The fact that the curve for the quadratically tapered layer 
lies below that of the homogeneous material confirms the improvement 
in performance, and we can now achieve a reflection coefficient of -50 
dB by choosing 0t/\ x ~ 0.64 corresponding to N = 9. Approximations 
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to t lif lone-dashed curve." in Figure > arc 

.it 

— = -u.oi lfij/id - o.or>i 

/\ j- 

A = 0.298 exp [3.263.*/ / A J 

where |/?| and A are as before. Compared with the homogeneous ma- 
terial the decrease in the number of elements required becomes more 
pronounced as |/?| is reduced. As previously mentioned, these- equa- 
tions hold for all real angles of incidence in the case of the anisotropic 
layer and for normal incidence if the layer is isotropic. 



4 THREE-DIMENSIONAL VERIFICA- 
TION 

As noted earlier, a PML is particularly attractive for terminating a fi- 
nite element mesh in the simulation of microwave circuits. For these 
applications a PML has an advantage over a traditional ABC because it 
does not require an a priori knowledge of the guided wave propagation 
constant. To demonstrate the applicability of the design equations in 
three dimensions, we consider a shielded microstrip line and a rectan- 
gular waveguide. 

The microstrip line has width w = 0.71428 cm, substrate thickness 
0.12 cm and relative permittivity c r = 3.2, and is enclosed in a metallic 
cavity whose dimensions are shown in Figure 9. It should be noted 
that the height of the cavity from the microstrip line is sufficiently 
large to suppress any reflections from the cavity walls. As a result, 
the characteristic impedance of the line should be that same as if the 
line was in free space. The microstrip line was terminated using a two- 
section homogeneous uniaxial absorber having material parameters ! r , 
Ji T in the upper section and 3.2f r , j! r in the lower section to match 
the substrate. The calculations were carried out at several frequencies 
using an FEM code [8] and we show the results for 4.0 GHz. At this 
frequency the element width was 0.05A o and a five layer absorber having 
a total thickness of t = 0.25A o was used. With a = 0 the computed 
reflection coefficient of the transmission line structure as a function of 
/? is shown in Figure 9. Recognizing that most of the field is confined 
to the substrate, we have X x = A 0 / v /q = 0.559A o . Using t = 1.789A* 
and N = 5 in the design formulas (9) and (10) gives (3 — 1.07 and 
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Figure 9: (a) Geometry for the microstrip line, (b) computed |/?| for 
the microstrip line at 4 GHz with a = 0. t = 0.25A o and A' = 5. The 
intersecting straight lines indicate the predicted values. 

|ft| = -41 dB (indicated on the plot by the vertical and horizontal lines, 
respectively). These values agree reasonably well with the numerical 
results shown in Figure 9 where the maximum absorption occurs for 
& ~ 1 and |/?| ~ -50 dB. The fact that the minimum |tf| is lower 
than predicted is, perhaps, not surprising. We recall that the design 
curves are based on the worst case, i.e. the polarization providing the 
largest minimum |/?|, and the curve in Figure 9 resembles more the H 
polarization curve in b igure 5 than the E polarization which constitutes 
the worst case. 

The other geometry considered is the rectangular waveguide shown 
in Figure 10. The elements are 0.5 cm bricks and the absorbing layer is 5 
cm thick (10 elements used) with a material parameter 6=1 —j/3. The 
reflection coefficient was computed at 4.0, 4.5 and 5 GHz for various 
values of (3 in order to determine the maximum absorption point. The 
resulting reflection coefficients are plotted in Figure 10. Using equations 
(9) and (10) with N = 10 yields /3t/\ x = 0.574 and \R\ = -50 dB. The 
vertical and the horizontal lines in Figure 10 indicate the location of 
these values. Once again, the agreement with the predicted values 
is good, with only a slight discrepancy in the value of /3t/\ x and a 
deviation of about 5 dB in the anticipated reflection coefficient. There 
are two points to keep in mind here. First, the design criteria have been 
applied to a non-normal incidence case in a three-dimensional setting. 
Secondly, the real part a of the material parameter b was set to 1, 
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Figure 10: (a) Geometry for the rectangular waveguide, (b) computed 
|/?| for the waveguide at 4,0 GHz ( ), 4.5 GHz (• ■ •) and 5.0 GHz 

( ) with q — 1, t = 5.0 cm and A = 10. The intersecting straight 

lines indicate the predicted values. 


demonstrating that the criteria may still apply when q ^ 0. 


5 CONCLUSION 

A uniaxial perfectly matched layer provides a powerful means for trun- 
cating finite element meshes close to the modeled structure. By prop- 
erly selecting the material properties and sampling rate, almost any 
desired level of absorption can be attained, and typically very few sam- 
ples (less than five) are needed to achieve a reflection coefficient of -40 
dB over a wide range of incidence angles. In this paper we described a 
detailed study of three types of layer material including homogeneous 
and inhomogeneous uniaxial ones, and by identifying the scalable pa- 
rameters of the layers, universal design curves and formulas were de- 
veloped. The curves or formulas can be used to specify the numerical, 
geometrical and electrical parameters of the PML to achieve a desired 
absorption down to -60 dB or lower. They are valid for all real angles of 
incidence for the anisotropic layer and restricted to normal incidence for 
the isotropic layer. As expected, a lower material loss requires a thicker 
absorber to produce the same reflection coefficient. On the other hand, 
a higher attenuation rate requires more samples to attain a lower re- 
flection coefficient in a numerical implementation. An inhomogeneous 
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(tapered | PML is l»et t or than a homogeneous sime ih<‘ material 

loss ran he increased to larger values close to the metal hacking where 
t he field is smallest . 

To test the applicability of the design criteria in a three-dimensional 
setting, a shielded microstrip line and a rectangular waveguide were 
considered. With both structures terminated with a homogeneous 
PML. the results were in reasonable agreement with prediction, and 
the discrepancies were no more than could be expected in view of the 

conditions under which the criteria were established. These conditions 
are: 


(i) use of a particular one-dimensional FEM code 

(n) based on the worst case polarization, i.e. the polarization for 
which the minimum reflection coefficient is largest 

(iii) assumption of a pure imaginary propagation constant, i.e, a = 0 , 
in the layer. 

Condition (iii) is a requirement for scalability, and though small 
values of q are still admissible, the condition is clearly inappropriate 
if there is substantial power at complex angles of incidence for which 
targe q is required for absorption. If the polarization can be specified, 
(11) is also inappropriate, and the design criteria may underestimate the 
performance that can be achieved. In any given problem where there 
is the luxury of testing a variety of layer specifications, it is probable 
that a performance can be achieved which is better than that predicted 
by the criterion, but even then the design values are a logical place to 
start. In the more likely situation where prior testing is not feasible, 
we believe that the design criteria provide a logical basis for specifying 
the parameters of the PML and its sampling. 
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Application and Design Guidelines of the 
PML Absorber for Finite Element 
Simulations of Microwave Packages 

J. Gong*. S. Legault*. A . Botros* and .I.L. Yolakis* 


Abstract 

The recently introduced perfectly matched layer(PML) uniaxial absorber for fre- 
quency domain finite element simulations has several advantages. In this paper we 
present the application of PML for microwave circuit simulations along with design 
guidelines to obtain a desired level of absorption. Different feeding techniques are also 
investigated for improved accuracy. 


I. Introduction 

In the numerical simulation of 3D microwave circuits using partial differential approaches, it 
is necessary to terminate the domain with some type of non-reflective boundary conditions. 
When using frequency domain PDE formulations, such as the finite element method, the 
standard approach is to employ some type of absorbing boundary conditions(ABCs) [1], [2], 
[3]. Also, the use of infinite elements [4] or port conditions [5] have been investigated. All of 
these mesh truncation methods require a priori knowledge of the dominant mode fields and, 
to a great extent, their success depends on the purity of the assumed mode expansion at 
the mesh truncation surface. Larger computational domains must therefore be used and the 
accuracy of the technique in computing the scattering parameters could be compromised. 

Recently, a new anisotropic (uniaxial) absorber [6] was introduced for truncating finite 
element meshes. This absorber is reflectionless(i.e. perfectly matched at its interface) for all 
incident waves, regardless of their incidence angle and propagation constants. As a result, 
it can be placed very close to the circuit discontinuity and is particularly attractive for 
terminating the computational domain of high density microwave circuits where complex 
field distributions could be present. 

Although the proposed uniaxial PML absorber has a perfectly matched interface, in 
practice a finite metal-backed (say) layer must be used which is no longer reflectionless due 
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to the presence of the pec (see Fie. 1 i. I, is t herefore of inters, n, opt imize . hr absorp. ivit v of 
J; M ; ■ - sH, ' r,lon of < parameters to achieve a civet, reflect ivity n il 1, a minimum 

layrr thickness. In this paper, we present guidelines fur implement ine the PML absorber to 
truncate finite element meshes in microwave circuit simulations. Lxample microwave circuit 
calculations are also given to demonstrate the accuraev of the PM 1. al, M »rl.er the FFM 
simulator and the feed model. More examples will be presented a. the conference. 

II. Absorber Design 

An extensive study was carried out using two-dimensional(see Fig. 1 ) and three dimensional 
models (see Fig. 2) in order to optimize the absorber's performance using the minimum 
thickness and discretization rate. As expected, the absorber's thickness, material properties 
and the discretization rate all play an equally important role on the performance of the P.ML. 
The t\pical field behavior interior to the absorber is shown in Fig. 3. As seen, for small A 
\alues the field decay is not sufficient to eliminate reflections from the metal backing. For 
arge 0 values, the rapid decay can no longer be accurately modeled by the FEM simulation 
and consequently the associated VSWR increases to unacceptable values. However, an op- 
timum value of 0 which minimizes the reflection coefficient for a given layer thickness and 
discretization can found. The parameters 0 and t play complimentary roles and the study 
shows that the PML absorber's performance can be characterized in terms of the product 
M (a scalable quantity when a = 0) and the discretization rate. A two-dimensional analysis 
was carried out to determine the optimum values of g and N (the number of samples in the 
PML layer) for maximum absorption near normal incidence. It was determined that given 
a desired reflection coefficient \R\ for the PML absorber, the optimum f and N values are 
approximately given by the expressions [7] 9 

0t 

= -0.01061/?! + 0.0433 




A r = 0.147 exp 


7.353^ 

*9 J 


where |/?| must be given in dB and N is equal to or exceeding the right hand value. As an 
example, if we desire to have a value of \R\ equal to -50 dB, from the above formulae we 
have that ^ % 0.58 and N = 10. It should be noted that though the design formulae were 
derived with a = 0 they also hold for small non-zero values of a. 


III. Feed Excitation 

Two feed models were used in conjunction with the scattering parameter extraction method. 
One was the horizontal current probes (Fig. 2) linking the back PEC wall with the beginning 
of a microstrip feed line. About 3 to 5 horizontal probes were needed for convergence and 
this scheme proved more accurate that the usual single vertical probe. 

The other feeding scheme employed here involved the specification of the quasi-static 
TEM mode at the microstrip line port. In the context of the FEM, the excitation is intro- 
duced by imposing boundary conditions across the entire cavity cross section through the 
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input port . Those conditions also serve to suppress backward relied ions from tiie modeled 
ciicuit discontinuity. C onsecjueni ly . they can l>e placc'd close to tin 1 discont mini v without 
compromising the accuracy of the scattering parameter extraction. 

IV. 3D Modeling Examples 

The PML performance' as predicted by t lie formulae was investigated bv using it to truncate 
the domain of 3-D microwave circuits. For example. Fig. -1 shows the optimum value of 
Xj zz 0.% obtained from the above design equations compares well with the results of the full 
wave FEM analysis of the microstrip line shown in Fig. 2. The 3-D FEM computat ions were 
carried out using A = 5 for modeling the PML absorber across its thickness and from the 
given formulae, it follows that R = — A\dB and this agrees well with the optimum value shown 
in Fig. 4. Another example is the meander line shown in Fig. 5. For the FEM simulation, 
the structure was placed in a rectangular cavity of size 5.8mm x IS. Omni x 3.175mm. The 
cavity was tessellated using 29 x 150 x 5 edges and only 150 edges were used along the y-axis. 
The domain was terminated with a 10 layer PML, each layer being of thickness t = 0.12mm. 
The 5’n results are shown in Fig 6 and are in good agreement with the measured data [8]. 
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Figure 1: Illustration of wave incidence upon a perfectly match interface (PML) with and 
without metal backing. 
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Figure 2: Shielded microstrip line terminated by a perfectly matched uniaxial absorber layer. 
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Figure 3: Illustration of the field decay pattern inside the PML layer. 



Figure 4: Reflection coefficient vs 2(3t/\ g with a=l, for the shielded microstrip line termi- 
nated by the perfectly matched uniaxial layer. 


95 



Figure 5 



Figure 6. Comparison of calculated and measured results for the meander line shown in 

Fig-5- 
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Abstract 

This chapter provides an overview of the finite element method (FEM) as applied to 
electromagnetic scattering and radiation problems. It begins with a review of the method- 
ology with particular emphasis on new developments over the past five years relating to 
feed modeling, parallelization and mesh truncation. New applications which illustrate the 
method’s capabilities, versatility and utility for general purpose applications are discussed. 
Specifically, new finite element modeling of antennas on doubly curved platforms, bandpass 
radomes and jet engine scattering are presented using a variety of mesh truncation schemes 
(boundary integral, absorbing boundary conditions and perfectly matched absorbers) are 
presented. Also, an entire section of the chapter is devoted to a reduced ordering method 
based on the Pade expansion. This reduce ordering method is used to obtain broadband 
responses from a few data points of the entire response. It is introduced in this chapter 
for the first time in connection with hybrid finite element systems and promises substantial 
computational savings. 


1. Introduction 

Over the past 10 years we have witnessed an increasing reliance on computational methods 
for the characterization of electromagnetic problems. Although traditional integral equation 
methods continue to be used for many applications, one can safely state that in recent years the 
greatest progress in computational electromagnetics has been in the development and application 
of partial differential equation(PDE) methods such as the finite difference-time domain and 
finite element (FEM) methods, including hybridizations of these with integral equation and high 
frequency techniques. The major reasons for the increasing reliance on PDE methods stem 
from their inherent geometrical adaptability, low O(N) memory demand and their capability 
to model heterogeneous (isotropic or anisotropic) geometries. These attributes are essential 
in developing general-purpose analysis/design codes for electromagnetic scattering, antennas. 
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.irrowavo circuits and biomedical applications. For example, modern aircraft geometries contain 
large non-metallic or composite material sections and antenna geometries mav involve many 
layers of materials, including complex microwave circuit feeding networks. Although, the moment 
method continues to remain the most accurate and efficient approach for sub wavelength hzc 
bodies of simple geometries. PDF methods and hybrid versions of these have shown a much 
greater promise for large scale simulations without placing restrictions on the geometrv and 
material composition of the structure. 

This chapter is a selected review of hybrid finite element methods and their applications to 
scattering and antenna analysis. We begin by first introducing the mathematical basics of the 
method without reference to a specific applications and in a manner that identifies the method's 
inherent advantages in handling arbitrary geometrical configurations which may incorporate 
impedance boundary conditions and anisotropic materials. In this section we also identify the 
key challenges associated with the implementation of the method such as mesh truncation and 
feed modeling for antenna applications. Section 3 of the chapter reviews the various mesh 
truncation schemes to be employed on the outer mesh surface for the unique solution of the 
vector wave equation. Absorbing boundary conditions, integral equations and artificial absorbers 
are discussed, all leading to different versions of hybrid FEM methodologies, and we comment on 
their accuracy and ease of implementation. Section 4 presents several approaches for antenna feed 
modeling in the context of the FEM, including coaxial as well as aperture coupled feeds. The next 
section is devoted to parallelization issues specific to finite element codes. We give performances 
of typical FEM codes and provide storage and implementation guidelines for maximizing code 
performance on parallel computing platforms. Section 6 describes a reduced order modeling 
approach for extrapolating broadband responses from a few data points. The method is based 
on the Asymptotic Evaluation Method(AWE) and is developed here for hybrid finite element 
systems. The final section of the paper gives some additional representative applications of the 
finite element method to scattering and antenna problems. 


2. Theory 

2.1. FEM Formulation 

Consider the antenna and scattering configurations shown in Figure 1. In the case of a scatterer, 
the entire computational domain is enclosed by a fictitious surface S a that may encompass a 
variety of composite/dielectric volumes as well as metallic, impedance and resistive surfaces. 
The antenna geometry is assumed to be recessed in some doubly curved surface. In this case, 
the bounding surface S 0 may either be located as shown or can be coincident with the antenna 
aperture. As usual, the recessed cavity is intended to house the radiating elements and their 
feeding structure such as coaxial cables, striplines, microstrip lines or aperture coupled feeds. 
The cavity may encompass any inhomogeneous or anisotropic material including resistive cards, 
lumped or distributed loads and so on. 

The goal with any finite element formulation is to obtain the solution of the vector wave 
equation 


V x 


/A 1 


(V x E) 


k 0 c r ’ E — jk 0 Z 0 J t 


v x (p ; 1 


■Mi) 


( 1 ) 
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Surface So enclosing 
the computational domain 





Finite 

element 

mesh 


Dielectric of 
volume Vd 
(enclosed by Sd) 


Resistive or 
impedance 
surface, S r 


(a) 


Conformal 




line Absorber 

termination 


(c) 


Figure 1: Computational domains for FEM analysis, (a) the various regions enclosed by So , (b) typical 
tetrahedral mesh, (c) computational regions for antenna analysis 
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m whirl, E reprints the total electric field. l?../7 r i deiiol r t he relative permittivity and perme- 
ahdity tensors of the computational domain. k b is the free '-pare wavenumber and i- the free 
spare intrinsic impedance. In addition. J, and M. denote ihe impressed elertrir and macetir 
sources, respectively, and represent the exc itation for the antenna problem. As is well known, 
the standard finite element (FE) solution scheme' is to consider the- weak form of ( 1 ). instead of 
solving it directly. This can be achieved by extremizing the pertinent functional [D 


F( E) = 


+ 

+ 

+ 


\ Jjf v {( V x E) • [K l ■ (V X E)] - k-(f T • E) e} dr 

JJf v E • UkoZoJ, + V X (jj;' . M .)] f/r 


jkoZo [ I E (H x i i ) ds 

JJSo+Sf 


\ ikA {SL> 


x E ) • ( n x E ) d, 


• + £JJL, 


EE dV 


where R denotes the resistivity or impedance of the surface 5 r . S f stands for the junction opening 
to possible guided feeding structures and H is the corresponding magnetic field on 5 0 and Sj 
whose outer normal is given by n. Also, V s is the source volume and VJ is the volume occupied 
by the load Z L , whose length and cross section are given by d and s, respectively. It is noted 
that Sr must be closed when it satisfies an opaque impedance boundary condition but can be 
open (i.e. a finite plate) if it satisfies the penetrable resistive sheet condition [2] 


h x n x E = -Rn x (H + - H - ) 

where H ± denote the fields on the two sides of the surface S r . As seen from (2), in spite of the 
discontinuity in the magnetic field, no special care is required for the evaluation of the integral 
over S T . The explicit knowledge of H is, however, required over the surface and Sf (referred 
to as mesh truncation surfaces) for the unique solution of E. 

As an alternative to (2), we can instead begin with (1) by weighting it with a testing function 
T. Subsequently, application of the divergence theorem yields the residual 

< R, T > = /// v {(VxT).[r.(VxE)]-a.E.T} rf F 

+ JJJ v T • [jk 0 Z 0 J, + V x (% l • M,-)} dV 

+ jk o Z 0 if T • (H x fi)dS 

V j So + S j 

+ 3hZ ° {Us. s ( " x T) ' ( " x E)ds + zls III T ■ E<n 1 (3) 

This is typically referred to as the weak form of the wave equation, whereas (2) is referred to as 
the variational form. It will be seen later that when set equal to zero, (3) yields the same system 
of equations as those deduced from the functional (2). Therefore both methods are equivalent. 
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7 hr funrtional ( 2 ) or the weighted residual ( 7 i can l>e cast into a discrete 1 svstem of rquat ions 
for the solution of E. 7o accomplish this, it is first necessary tu subdivide tin* volume' as a 
collection of small elements such as those shown in Figure 2 [:*j. Within each volume element, 
the field can then he expanded as 


Right Angled Skewed Brick Curvilinear Brick 

Brick 


Tetrahedron Distorted Prism Cylindrical Shell 

Figure 2: Illustration of the various elements for tessellating three dimensional volumes 


M=# of edges 

E‘= ^ E{Wl = [W'f {£'} (4) 

k= 1 

in which W£ are referred to as the edge-based shape or basis functions of the eth element in the 
computational domain. In contrast to the traditional node-based shape functions, the edge-based 
shape functions are better suited [5] for simulating three dimensional electromagnetic fields at 
corners and edges. Moreover, edge-based shape functions overcome difficulties associated with 
spurious resonances [6]. They were proposed by Whitney [7] over 35 years ago and were revived 
in the 1980s [8], [9]. Their representation is different for each element but have the common 
properties [10] of being divergence free (i.e., V • W fc = 0 within the element) and normalized in 
such a manner so that the expansion coefficients El represent the average field value across the 
kth edge of the eth element. One disadvantage of the edge-based elements is that they increase 
the unknown count. However, this is balanced by the increased sparsity of the resulting stiffness 
matrix. A detailed mathematical presentation of the edge-based shape functions for various two 
and three dimensional elements (bricks, hexahedra and tetrahedra) [11]-[14] is given by Graglia 
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(i 


(I. al. [15] in this issue of the Transactions. Linear as well a- higher order (curvilinear 
are discussed [lGj-[l9j. 


element s 


2.2. Discretization and System Assembly 

Once the computational domain has been tessellated using appropriate elements and shape func- 

•v, 

tions. the discretization of (2) or (3) proceeds by introducing the expansion E = ^E'. with 

E as given b\ (4) and A t denoting the number of elements comprising the computational do- 
main. For the functional (2). the system of equations is constructed by setting (Rayleigh-Ritz 
procedure) 


dF{ E) 

dEi 


= 0 , 


f = 1.2 .Y e ; k =1.2 M 


(5) 


whereas in the case of (3) Galerkin’s approach is used by setting T = W£. Regardless of the 
procedure, the resulting system will be of the form 


EmtC'} + £>*]{£•} + £>"} + ^{C-} = 0 (6) 

e=l 5=1 e= 1 5=1 

in which the brackets denote square matrices and the curly braces refer to columns. Among 
the various new parameters, [ A e ] is referred to as the element matrix and results from the 
discretization of the first volume integral in (2) or (3); N s is equal to the aggregate of the 
surface elements (quadrilaterals for hexahedra or triangles for tetrahedra and prisms) used for 
the tessellation of S 0 , Sj and 5 r ; the column {/\ e } results from the discretization of the source 
integral over V, ; [fi a ] is the matrix associated with the surface integrals over So + 5/ in (2) or 
(3) with TV,, aggregated surface elements; and finally {C*} results from the discretization of the 
corresponding surface integral in (2) or (3) involving the external incident field. Basically, {C 5 } 
provides the excitations for scattering problems whereas {A e } is the corresponding excitation 
column for internal sources as is the case with antenna problems. 

The entries of the element matrix [A e ] are specific to the choice of the shape functions and 
are compactly given by 

Ml = Jjf v {(£W'P r ]-’[£WT - WW'^WW'] 1 } dV (7) 

with 


[DW e ] T 


iWY-i-AK} 

1 


( 8 ) 


and V e denotes the volume of the discrete element while the subscripts (x,y,z) in (8) imply the 
{x,y,z) components of the shape function W e . If lumped loads are present (i.e., in the presence 
of a volume integral over V^), the diagonal entries of [zl e ] are simply modified with the addition of 
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tli( torn) jk.,Z,-.d /Zj . I erfectb conducting posts- located at tlio /.■' edge arc handled lu setting 
tho corresponding total fields equal to zero and rearranging the final matrix as discussed later. 
I he internal source excitation column is simply given by 



The specification of the matrix [ B s ] and the entries of the corresponding excitation column 
{C 5 } due to the incident field can not be completed without first introducing some a priori 
relationship between the E and H fields on .S 0 . This relationship (or boundary condition on 
S 0 ) is referred to as the mesh termination condition and its form depends on the physics of 
the problem. For example, in the case of problems where the entire computational domain is 
enclosed by a perfect electric or magnetic conductor, the unknown column {£*} is eliminated 
from {E}. For scattering and radiation problems, the mesh termination condition should be 
such that the artificial boundary S D appears transparent to all waves incident from the interior 
of the computational domain. Clearly, if S 0 is placed far enough from the scatterer or radiator, 
the simple Sommerfeld radiation condition provides the appropriate relationship between E and 
H. However, this is not practical since it will yield large computational domains. To bring S 0 as 
close as possible to the scattering or radiating surface, more sophisticated boundary conditions 
must be introduced on S 0 . These mesh termination conditions are critical to the accuracy and 
efficiency of the formulation and are some of the major bottlenecks in the implementation of the 
FEM. Regardless of the employed mesh termination approach, after carrying out the sums in 
(6), a process referred to as matrix assembly, the resulting matrix system will be of the general 
form 



In this, {E v } denotes the field unknowns within the volume enclosed by S Q + Sj whereas { E B } 
represents the corresponding unknowns on the boundaries S a and/or Sj. The excitation column 
{6 V } results from the assembly of {A' e } and similarly {6 B } is associated with {C*}. The matrix 
[.4] is very sparse (9 to 30 non-zero entries per row) and this is a major characteristic and 
advantage of FEM. By using special storage schemes and solvers suitable for sparse systems, 
the CPU and memory requirements are maintained at O(N) and scalability can be attained on 
multiprocessor platforms. 

The boundary matrix [Q] in (10) is associated with the boundary fields {E B } and its specific 
form is determined by the employed mesh termination condition on S 0 . Over the past five years, 
much research on FEM has concentrated on the development of mesh truncation schemes which 
minimize the computational burden without compromising accuracy. As will be seen later, this 
compromise is difficult to attain and is subject to the computational problem being considered. 
For the purpose of our discussion, we will subdivide the various mesh termination schemes into 
three categories: (1) exact boundary conditions, (2) absorbing boundary conditions (or ABCs), 
and (3) artificial absorbers. Exact boundary conditions provide an integral relation between the 
electric and magnetic fields, and because they are exact, they permit the placement of S a very 
near or exactly on the surface of the scatterer or radiator. The resulting formulation is referred 
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to as the finite element-boundary integral (FE-BI i method and combines the adaptability of the 
FEM and the rigor of the boundary integral methods. However, it yields a fully populated matrix 
[Q\ and as a result, the FE-BI method is associated with higher computational demands even 
though only a small part of the overall system is full (see Figure :1a). To alleviate the higher 
computational demands of the rigorous FE-BI method. ABCs or artificial absorbers (AAsfcan 
instead be used to terminate the mesh. Both of these are approximate mesh termination methods, 
but lead to completely sparse and scalable systems of equations. Basically, the sub-matrix [£] 
is eliminated with the application of ABCs or AAs (see Figure 3b). In the case of ABCs, a 
local boundary condition in the form of a differential equation is applied on to relate E and 
H so that S c appears as transparent as possible to the incident fields from the interior. The 
resulting method will be referred to as the finite element-ABC (FE-ABC) method and has so 
far been the primary method for general-purpose scattering computations. Finally, the use of 
artificial absorbers (including perfectly matched layers or PMLs) [20] have recently gained major 
attention because of their potential for greater accuracy and inherent implementation simplicity. 
In the context of the finite element-artificial absorber (FE-AA) method, the mesh is terminated 
by using a material absorber (typically non-realizable in practice) to absorb the outgoing waves 
and suppress backward reflections. Below, we briefly discuss the specifics for each of the three 
mesh termination schemes. 


3. Mesh Termination Schemes 

3.1. Finite Element-Boundary Integral Method 

The FE-BI method was introduced in the early seventies [21], [22] as a natural extension of 
the FEM for modeling unbounded problems. However, because of its larger computational 
requirements, the method did not enjoy a widespread application to electromagnetics until the 
late eighties [23, 24]. In accordance with the FE-BI method, the relation between E and H on 
So is provided by the Stratton-Chu integral equation 

n x H(r) = h x H' nc (r) 

+ fix If {[h' x H(r)] x V'G 0 (r, r') + jk 0 Y 0 [h' x E(r)]G c (r, r') 

J J S 0 

+ j^V.|»'xE(r')]VC,(r,t')J iS' ( 11 ) 

where r and r' denote the observation and integration points, respectively, and G 0 (r, r') = 
ex p(~jk 0 | r — r |/ (47r |r — r ; |) is the free-space Green’s function. The above is the most general 
form of the boundary integral and places no restrictions on the shape of S 0 with the exception 
that it must be closed as shown in Figure 1(a). Provided S 0 is placed just above the outer 
boundary of the scatterer or radiator, any material composition which is interior to S a can be 
handled with ease using the FEM without relation to the boundary integral. This form of the 
FE-BI was used by Yuan [25] and later by Angelini et. al. [26], Antilla and Alexopoulos [27], 
and Soudais [28] to model three dimensional scatterers with anisotropic treatments. The method 
has also been successfully used for biomedical simulations [29], [30]. 
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Id 


The discrelization of 


1 1 j follows by introducing the expansion unci the s Ih surface element 


r. = # of 

E(r)= V E^(R) = [Nf {/;*} ,!■_>! 

*T=l 

with a similar expansion for the H field (if necessary). The surface shape functions S‘(r) can be 
set equal to W e (r) when the latter are evaluated on the surface of the volume element . For the 
tetrahedron, n x S* are simply equal to the traditional basis functions given by Rao ft. al. [31], 
However, as pointed out in [32], S s can be chosen independently of W' provided care is exercised 
when (11) is substituted/combined with (2) or (3). Regardless, when (11) is substituted into 
(2) or (3) after discretization we will get the partly sparse, partly full system [33] given in (10) 
and illustrated in Figure 3(a). Finally, we point out that since S 0 is a closed surface, the final 
system is not immune to false interior resonances due to the boundary integral. In this case, the 
combined field formulation [34] or the simpler complexification scheme [35] may be employed. 

When S 0 coincides with an aperture in a ground plane (see Figure 1(b)), the integral equation 
(11) can be simplified substantially [36]. Specifically the integral relation between E and H on 
the aperture reduces to 


n x H = n x H 5 ° + 



h x G(r, r') [E(r) x n] dS 


(13) 


where G is the dyadic Green’s function of the second kind, with h x V x G = 0 on S c and for 
planar platforms it reduces to , 



g-jMr-r'l 

47r(r — r'| 


(14) 


In this, I is the unit dyad and the factor of 2 is due to image theory. Also, H s ° is equal to the 
sum of the incident and reflected fields (in the absence of the aperture) for scattering or zero for 
antenna parameter computations. For non-planar S 0 , HU° is equal to the field scattered by an 
otherwise smooth surface again with the aperture removed. In this case, the Green’s function 
should also be modified accordingly with respect to the non-planar platform. One of the first 
implementations of the FE-BI for radiation and scattering from rectangular apertures/antenna 
recessed in a ground plane was given by Jin and Volakis [14], [33], [36] and was later extended to 
antenna analysis on planar [37], [38], [39] and cylindrical [40] platforms. 

The FE-BI method is particularly attractive in terms of CPU and memory requirements 
when the [£] matrix is Toeplitz and can therefore be cast in circulant form [41], [42], [43], [44], 
In this case, the entire system can be solved using an iterative solver [45] in conjunction with 
the FFT to reduce the CPU requirements down to 0(N se log N se ) for the boundary integral 
sub-matrix. The FFT is simply used to carry out the matrix-vector product [Q] { E B } within 
the iterative solver. For example, if the symmetric biconjugate gradient (BCG) method [45] 
is employed to solve (10) and rectangular elements are used, the storage requirement is only 
4Af ee + 8 N se , where N ee and N 3e denote the number of edges within the computational volume 
and on S a , respectively. For triangular surface elements, the storage requirement is about 4.5 
times larger due to the presence of the diagonal edges across the quadrilaterals. Whether the 
full matrix [Q] is Toeplitz or not depends on the shape of the surface S a and on the uniformity of 
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tlic mesh. It can Ik- shown that [Cj will he Toeplitz for planar and cylimlrical surfaces provided 
the surface grid is uniform. The Toeplitz form of [C] has been exploited with great success ami 
systems with more than 0.5 million unknowns spanning apertures of 15A • 15A have been solved 
accurately on a desktop workstation. For example, we observed that a ISO. 000 unknowns system 
converged in 57 iterations and for another system of 125.000 unknowns convergence was achieved 
after 60 iterations with an average CPF time of 2 sec/iteration on an 1IP0000/75O rated at 
24.7 M Flops. Because of these impressive performances, it is advantageous to transform [i7] to 
a Toeplitz matrix even when the surface grid is not uniform. To do so. a uniform grid can be 
superimposed onto the non-uniform mesh (see Figure -1) and the edge fields on the two grids can 
be related via a connectivity matrix [39]. In this manner, non-rectangular antenna elements and 
apertures can be modeled with the full geometrical adaptability of the finite element method and 
without compromise in accuracy. It is important to note that similar transformation matrices 
can be exploited for decomposing computational domains as done, for example, in [97]. 



Figure 4: Overlay of a uniform grid on an unstructured mesh for implementation of the FFT to perform 
the matrix-vector products. 

When So is not planar, the boundary integral matrix will inevitably cause large CPU and 
storage requirements for large N 3e to the point where F E-ABC or FE-AA methods become 
more attractive at the expense of accuracy. Recently, though, techniques have been proposed 
which show promise in reducing the computational requirements of the matrix-vector product 
\Q\ {F }. Among them, the fast multipole method [48], [49] can reduce the operation count in 
carrying out the matrix-vector product from O (N* e ) down to O (N$ e 5 ) and reductions down to 
0{N} e 33 ) have also been proposed. The main idea of the FMM [49] is to subdivide the surface 
So into groups, each containing approximately M ae « y/NZ unknowns. By rewriting the free 
space Green’s function as an expansion [50] or by introducing its far field approximation [51], 
it can be shown that the interactions between the unknowns within the groups can be carried 
out in O ( N se M se } operations whereas the interactions between groups can be carried out in 
0(N? e /M se ) operations. The sum of these two operation counts yields a total operation count 
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of ()( A' 3 V’) for the Foundary integral matrix when \l s , ^ x \„. and this must he a tided u. the 
(9(A'„ ) operation count associated with the sparse matrix [A]. Details for the implementation 
of the FMM are given in [4S]. [-19] and when combined with the FEM. we can refer to it a- 
the FE-FMM method [52]. A comparison of the FE-BI (with LE and conjugate gradient- EFT 
solvers) and the FE-FMM methods for the calculation of the scattering by a large groove in a 
ground plane is tabulated in Table E The groove was 50A long and 0.35A deep (all metallic) 


Groove Width 

Total Unknowns 

BI Unknowns 

CPI Tunt for HI (Minute s. seconds ) 

FE-BI (CG) 

FE-FMAM 

FE-BI (C'GFPT) 

25 A 

2631 

375 

(8.48) 

(3.26) 

(1.41) 

35 A 

3681 

525 

(16.34) 

( 5.55) 

(3.24) 

50A } 5256 

750 

(45,1) 

(14.31) 

(5.40) 


Groove Width 

Storage of BI (KB) | 

FE-BI (CG) 

FE-FMM 

25 A 

1072 

277 

35A 

2102 

458 

50 A 

4291 

T 

oo 



B MS error (dB ) 

Groove Width 

Fe-fmm 

25 A 

M2 

35 A 

1.2 

50A 

1.36 


Table 1: CPU Times, storage requirement and error 

and filled with a dielectric having c T = 4 and ^ i r = 1. By using a sampling rate of about 15 
edges per wavelength, the system unknowns were N ee = 5256 and N 5e = 750. From the data 
in Table 1, it is clear that the FE-FMM is more than 3 times faster and requires 3-5 times less 
memory than the FE-BI when LU decomposition is used for solving the FE-BI system. Because 
of the grouping/averaging of the surface elements, the FE-FMM exhibits certain errors which 
are functions of many parameters [53] and these errors must be kept in mind as the system size 
is increased. It is expected that the FMM will exhibit greater speed-ups as the system size is 
increased. Nevertheless, the special implementation of the FE-BI using the FFT is still by far 
the most accurate approach. 

We close this section by noting that another approach for treating the boundary integral 
matrix- vector products is via the adaptive integral method (AIM) [54], This approach relies 
on the introduction of multipole expansions (similar to Taylor series expansions) to replace the 
sources on S 0 by equivalent point sources placed on rectangular grids. In this manner the CPU 
requirements can be reduced by making use of three dimensional FFT routines. It can again 
be shown that AIM reduces the boundary matrix CPU requirements down to O (N 15 ) or even 
down to 0 ( A^ 33 ). 

3.2. Finite Element-Absorbing Boundary Condition Method 

The goal with any ABC is to eliminate backward reflections from S 0 , and a variety of these 
boundary conditions have been proposed over the years beginning with those of Bayliss and 
Turkel [55], and Engquist and Majda [56]. More recently, other ABCs such as those by Webb 
and Kanellopoulos [57], and Chatterjee and Volakis [4] (see also Senior and Volakis [2]) have 
been proposed. All ABCs provide an approximate relation between the E and H fields on the 
surface S 0 which in most cases is derived by assuming a field expansion in inverse powers of r, 
where r is the radial distance from the center of S 0 . If the ABC annihilates the first (2m + 1) 
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inverse powers of r. it is then referred to as an mth order ABC. The second order ABC derived 
by Ohatterjee and \olakis [-4] is given by 


where 


jkuZ 0 (W 


h) = a ■ E + j ■ X > f,(Y >. E Mt ) u \ - -• • V,i V EN ;I T 


ii-Vi 


0 = 
$ = 

1 = 


Ql 
O 2 

Ak 

Z) 


Ql / 1 / 1 + 02^2^2 

^(Ml + ^ 2^2 ) 

• i ( J *0 “I" 3K m — — 2 K] ) / 1 / 1 H — — (jA.‘o + 3 at T 

J™ 0 ^ m J KQ 

T )[‘^ K m ~ K 9 + D(jk 0 ~ a T ) + K j + K m Ak] 

7 )[^ K "m — K S T D(jk 0 — K 2 ) + Kt — K m A«] 

1 


— 2 k_>)GC 


D - 2ac, 

AC ] — K 2 


— 2 7T0 -F 5 k, 


In the above equation, the subscripts n and f denote surface normal and tangential components 
of the field, / 1 , ^ 2 a r^ the principal surface directions and Ki,K 2 are the corresponding curvatures, 
K m = («i + ^ 2)/2 is the mean curvature and K g = k\k 2 is the Gaussian curvature. For Kj = k 2 
(spherical termination boundary), (15) reduces to the ABC derived by Webb and Kanellopoulos 
[57]. The latter authors have recently proposed a correction [58] for the implementation of these 
ABCs which should be incorporated in existing codes. 

Upon substituting (15) into the surface integrals in (2) or 3), separating the scattered field 
from incident field, and making use of vector differential and integral identities [ 59 ], we have 


jk 0 Z 0 JJ s E scot • (H x h)dS = JJs M£ ( 7 ‘) 2 + a 2 (E“ at ) 2 )dS 

+ // „(V x E scat ) 2 n dS 
- JJ s (V • E* C “‘)[V • (7 • E*“‘) t ]d5. 


(16) 


Next, expansion (4) or ( 12 ) is introduced and (16) is differentiated as in (5) to obtain the finite 
element equations ( 6 ) . Since the ABC is a local condition, the sparsity of the finite element 
matrix is preserved and the resulting system is again given by (10) with [Q] removed. However, 
in the case of arbitrarily curved boundaries, the matrix will not be symmetric except for planar 
and spherical boundaries. The system will be also symmetric for cylindrical boundaries only if 
linear edge-based expansion functions are used. These ABCs have been extensively validated 
for a number of antenna and scattering configurations^], [60]. Higher order ABCs have been 
proposed for a more effective suppression of the outgoing waves. These are often difficult to 
implement but as demonstrated by Senior et. al. [61], they provide greater accuracy and effort 
being placed so much closer to the scatterer or radiator. 
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Figure 5: Illustration of an antenna terminated by an isotropic absorbing material 


3.3. Finite Element-Isotropic Absorber Method 

In accordance with this method, the outgoing waves are suppressed by an absorbing dielectric 
layer placed at some distance from the antenna as shown in Figure 5. Typically, the layer is 
backed by metal and the finite element region extends all the way to the metal at which point 
the mesh is terminated by setting the tangential component of the total electric field to zero. 
This is equivalent to removing the integral over S 0 in (2) or (3). Obviously, the accuracy of 
the technique depends on how well the metal-backed absorber suppresses the incoming waves 
without introducing backward reflections. The plane wave reflection coefficient of the absorber 
shown in Figure 5 has been minimized over the entire visible angular range [46]. The fact 
that the dielectric has the same relative permittivity and permeability ensures that there is no 
reflection from the air-dielectric interface at normal incidence (perfect impedance match), but 
the performance degrades away from normal with the reflection coefficient reaching unity at 
grazing. This limitation led researchers to consider perfectly matched interfaces as discussed 
next. Nevertheless, even though the isotropic absorbers are not perfect, they are still useful in 
modeling antennas on doubly curved platforms. Figure 6(a) displays a sectoral microstrip patch 


Metallic 

Ground 

Plane 


Patch 




(a) (b) 

Figure 6: Sectoral microstrip patch on cone 

printed on a conical surface and Figure 6(b) shows the antenna resonance behavior for different 
cone angles. We clearly observe that the resonance frequency drops with the cone angle for the 
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same patch dimensions. It should he also remarked that the computed resonance frequency is 
within 3.2 percent of that predicted by the approximate cavity model, and this is reasonable. 

In generating the results given in Figure (>. the computational domain was discretized using 
2 358 Prisms resulting in 3.790 degrees of freedom. One frequency run took o.*> minutes on an 
HP9000 (Model 713/64 ) workstation rated at 22 MFlops. 


3.4. Finite Element- Anisotropic Absorber Method 

This mesh termination method is similar to that described in the previous section except that 
the layer is comprised of an anisotropic lossy dielectric which has zero reflection coefficient at 
the air-dielectric interface over all incident angles and can therefore be considered as a perfect Iv 
matched layer (PML) [20]. The intrinsic parameters of the PML are 


_ = 
a 


where a - ( Z/Z 0 ) 2 , Z is the intrinsic impedance of the medium being terminated by the absorber, 
whereas a, fi and the thickness of the layer are parameters which can be optimized for maximum 
absorption. 

The advantage of the anisotropic (over the isotropic) artificial absorber in terminating finite 
element meshes is illustrated in Figure 7 [62] where it is shown that the anisotropic absorber 
retains its low' reflectivity at oblique incidences except near grazing (</> = 0°). The performance 
shown in Figure 7 is based on a purely theoretical analysis but as shown in Figure 8, this 
performance can be realized in numerical simulations with careful choices of fit and the sampling 
rate N within the absorber. It has been found [63] that the reflectivity curve shown in Figure 8 is 
typical to most situations and the location of the minimum reflection coefficient can be predicted 
a priori. An extensive numerical study based on a two dimensional planar absorber model 
demonstrated that fit/ X x (A* = X 0 / cosep) is a scalable parameter and that given the desired 
reflection coefficient |/?|, the formulas [62] 


Q — j3 0 0 

0 q — jfi 0 

0 0 — F 


— = -0.0106|fl| +0.0433 

A -r 


N = 0.147 exp 


7 



(U) 

(18) 


can be used to choose fit/X x and the minimum number of discrete samples to achieve the desired 
absorption. In these, \R\ is specified in dB and A 0 denotes the free space wavelength. These 
formulas are in agreement with the actual numerical results in Figure 8 but it has yet to be 
determined how well they apply for curved perfectly matched layers which are placed conformal 
to scattering and radiating surfaces. Improvements to their absorptivity though can be attained 
by considering tapered layers and formulas similar to (17)-(18) are given by Legault [62] for one 
such tapered anisotropic absorber. 
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Figure 7: Reflectivity of a A/4 thick planar metal-backed perfectly matched absorber (/? = 1 — j2) as 
a function of incidence angle, (a) Geometry, (b) plane wave reflection coefficient vs. angle 


-25 
-30 
| -35 
g-40 
-45 

absorber 

1 * 3.09S cm * Microstrip -50 

0.75 LOO 1.25 1.50 1.75 2.00 

P 

(a) (b) 

Figure 8. Reflection coefficient of the PML for terminating a microstrip line as extracted from a 
numerical implementation 
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4. Feed Modeling 

For .scattering problems where the plane wave incidence is usually the 'source'', the richt hand 
side excitation has been explicitly given in (10) and will not be discussed further. However, 
for antenna parameter computations, the explicit form of { A ' } in (ft) will depend on the type 
of feeding scheme being employed. Below we discuss specific forms of {K' } corresponding to 
different feeding choices. 


4.1. Simple Probe Feed 

For thin substrates the coaxial cable feed may be simplified as a thin current filament of length 
1 carrying an electric current / /. Since this Filament is located inside the cavity, the first term 
of the integral in (2) or (3) needs to be considered for this model. Specifically, the ith (global) 
entry of the excitation vector A, becomes 

A, = jk 0 Z 0 I 1 ■ W,(r), i - ji.j 2 j m ( 19 ) 

where r is the location of the filament, m is the number of (non-metallic) element edges and j m 
is the global edge numbering index. In general, m such entries are associated with m element 
edges, and thus the index i goes from j x up to j m . This expression can be further reduced to 
A, = jk 0 Z 0 I /, provided that the tth edge is coincident with the current filament. 


4.2. Voltage Gap Feed 

This excitation is also referred to as a gap generator and amounts to specifying a priori the 
electric voltage V across the opening of the coax cable or any other gap. Since V = E • d, where 
J is a vector whose magnitude is the gap width, and E the electric field across the gap, we have 

that E{ — ^ cos q ’ w here cos0, is equal to 1 if the zth edge is parallel to d. Numerically, this gap 

voltage model can be realized by first setting the diagonal term A it of the {A} matrix equal to 
unity and the off-diagonal terms A tj ( i / j) to zero. For the right-hand-side vector, only the 
entry corresponding to the ith (global) edge across the gap is specified and set equal to the value 
Ei whereas all other entries associated with edges not in the gap are set to zero. 

4.3. Coaxial Cable Feed Model 

The simple probe feed of the coaxial cable is accurate only if the substrate is very thin. For 
thicker substrates, an improved feed model is necessary and this can be achieved by evaluating 
the functional 

F c = -jk 0 Z 0 J j (E x H) • zdS (20) 

over the aperture Sj of the coax cable. Assuming a TEM mode across 5/, the fields within the 
cable may be expressed as (see Figure 9) [64] 
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with 



I - 1 I 


h 


a — 



In these expressions. t TC is the relative permittivity inside the cable. E and H are the electric 


cavity 


patch 



cavity-cable junction 
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k 
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Figure 9: (a) Side view of a cavity-backed antenna with a coax cable feed; (b) Illustration of the FEM 
mesh at the cavity-cable junction (the field is set to zero at the center conductor surface). 


and magnetic fields, respectively, measured at z — 0 and Iq is the center conductor current. Also, 
{r,<t>,z) are the polar coordinates of a point in the cable with the center at r = 0. We observe 
that (22) is the desired constraint at the cable junction in terms of the new quantities ho and eo 
which can be used as new unknowns in place of the fields E and H. 

However, before introducing F c into the system, it is necessary to relate e 0 and h 0 to the 
constant edge fields associated with the elements in the cavity region which border the cable 
aperture. Since the actual field has a 1/r behavior in the cable, we find that 


AV = Ei(b-a) = e 0 /n-, i = N p {p = 1,2,..., N c ) 


(23) 


where AV denotes the potential difference between the inner and outer surface of the cable and 
N p denotes the global number for the edge across the coax cable. When this condition is used in 
the functional F c , it introduces the excitation into the finite element system without a need to 

extend the mesh inside the cable or to employ a fictitious current probe. The derivation of — — 
. . . dE i 

and its incorporation into the system is then a straightforward task [64]. As can be expected, 

the above feed model assumes the presence of only the dominant(TEM) mode at the cavity-cable 
junction, an assumption which may not be suitable for certain applications. Of course, the model 
can be improved by extending the mesh (say, a distance d) into the cable. The equi— potential 
condition will then be applied at z=-d, where all higher order modes vanish. 
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4.4. Other Feed Models 

liierc arc a few other commonly used feed models for simulating antennas and the associated 
network in the context of the finite element methods. In certain cases, tlx- structures mav contain 
detailed geometries which must be modeled with care to ensure the efficiency and accuracv of 
the simulation results. For instance, the configuration of an aperture coupled microstrip antenna 
may be efficiently modeled by applying the equi-potential continuity condition and the interested 
readers are referred to [G5] for details. Also in modeling microwave circuits as antenna fet'd 
network, the excitation location along the network may have to be placed far from the antenna 
for probe models, and thus the modal excitation is an alternative to the probes described in 
section 4.1 and 4.2. This reduces the size of the computational domain without compromising 
accuracy. The modal field distribution is typically obtained using a simplified analysis model to 
truncate the 3D FEM domain. A 2D FEM code can be used as well for geometries having the 
same cross-section as the original feed network. 

In general, the antenna feed or feed network can be accurately modeled in the context of the 
FEM. Moreover, unlike the method of moments (MoM), the FEM provides the field distribution 
in the entire 3D computational space and this is particularly useful for visualization around the 
feed region and on the antenna. 


5. Parallelization 

When considering 3D problems of practical interest, the unknown count of the computational 
domain can easily reach several million degrees of freedom. The sparsity of the FEM system 
(particularly for the FE-ABC and FE-AA methods) makes possible the storage of such large scale 
problems but even at 0(N) computational demands, their practical solution requires efficient use 
of parallel and vector platforms. Modern computing platforms can now deliver sustained speeds 
of several GFlops and CPU speeds in the Tflops range are within our reach. The inherent 
sparse matrices of PDE methods are particularly suited for execution on multiprocessor and 
vector platforms but the exploitation of these processors requires special storage schemes and 
techniques to perform the matrix-vector product required in the iterative algorithms at the Flop 
rates sustained on these multiprocessors. 

To parallelize and vectorize the FEM codes, it is essential to first optimize the execution of 
the iterative solvers which typically take-up 90% of the CPU time. Among them, the conjugate 
gradient algorithms (CG, BCG, CGS and QMR) have been found very attractive and a brief 
comparison of the pros and cons for these is given in [45]. The Generalized Minimal Residual 
Method (GMRES) is another iterative solver which can exhibit faster convergence rates. How- 
ever, it stores the direction vectors and as a result it requires much higher storage. For the 
discussion below we will primarily concentrate on the BCG and QMR algorithms and we note 
that the symmetric form of BCG requires minimal number of arithmetic operations (see Table 
2). A disadvantage of the BCG is its erratic convergence pattern whereas the QMR has smooth 
and monotonic convergence. However, neither BCG nor QMR can guarantee convergence and 
typically they both converge or not for the same problem. When considering the parallelization 
of a completely sparse system such as that resulting from the FE-ABC method, the following 
issues must be addressed: 
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( 'uin plex c 


Operat ion 


f | 

Matrix-vector Products 

nze 

nze-N 

! Vector Updates 

-IN 

3N 

Dot Products 

3N 

3X | 

Total # of Operations 

29 X 

3N 


N = # of unknowns 

nze = # of nonzero matrix elements 


Table 2: Floating Point Operations of BCG Per Iteration 


5.1. Storage of Sparse System 

The performance of the code is strongly dependent on the employed storage scheme. Since a 
typical FEM matrix has about 8.5 A« or so non-zero entries, it is essential that the non-zero 
elements be stored in a manner that keeps the storage requirements nearly equal to the non-zero 
entries and minimizes inter-processor communications. The ITPACK [66] and the compressed 
row storage (CRS) schemes are most appropriate for parallel computing. The ITPACK storage 
format is most efficient for retrieving the matrix elements and certainly the choice method when 
the number of non- zero elements are nearly equal for every matrix row. Basically, the ITPACK 
format casts the FEM matrix in a smaller rectangular matrix having the same rows as the original 
matrix and can be unpacked by referring to a pointer integer matrix of the same size. However, 
this rectangular matrix can contain as much as 50% zeros which results in space wastage. By 
using a modified ITPACK scheme, space wastage can be reduced down to 30%. Even with less 
wastage, the CRS format may be the most efficient storage scheme with some compromise in 
CPU speed. It amounts to storing [.4] as a single long row which can be uncompressed using two 
integer pointer arrays. For the symmetric BCG algorithm, the CRS format results in only 8.5 
N complex numbers and N integers. However, it should be pointed out that the CRS format is 
not appropriate for vector processors such as the C-90. For vectorization, it is best to organize 
the storage in sections of long vectors and to achieve this for our type of matrices the jagged 
diagonal format [67] appears to work best. Using this format the rows are reordered so that the 
row with the maximum number of non- zeros is placed at the top of the matrix and rows with the 
least non-zero entries are shuffled to the bottom. This reordering greatly enhances vectorization 
because it allows tiling of the shorter rows to yield very long vector lengths in the matrix- vector 
multiplication phase. Specifically, for some problem the jagged diagonal storage format allowed 
the matrix- vector multiplication routine to run at about 275 MFlops on a Cray C-90 whereas the 
same routine ran at 60 MFlops using the CRS format. The dot product speeds and the vector 
updates reached 550 MFlops and 600 MFlops for the same problem. Table 3 provides a relative 
comparison of CPU estimates on various computers. 

5.2. Interprocessor Communications 

For distributed memory platforms, the method of partitioning the stiffness matrix [.4] among 
the processors, the chosen storage scheme and the inherent unstructured sparsity of [.4] are all 
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Processors 
4 

1 £ of Processors. P 

; I iom i //-seo • iteration ' unknown i j 

1 Cray C-90 

1 ( 277 M Flops) 

(T T 

| KSR 1 

28 

1 .28 


58 

0.57 

Intel Paragon 

8 

3.42 

16 

1 .99 


32 1 

1.38 

[ IBM SP-1 

4 

1.47 


Table 3: CPF Time Per Unknown for Solving Typical FE-ABC Systems 


crucial to the overall speed of the code. An approach that has worked well on massively parallel 
processors (such as the SP-2. Intel Paragon, Convex Exemplar) is that of assigning each processor 
a section of the matrix and by dividing the vectors among the P processors. Thus, each processor 
is responsible for carrying out the matrix-vector product for the block of the matrix it owns. 
However, the iterate vector is subdivided among all processors, and therefore narrow-band or 
structured sparse matrices have an advantage because they reduce interprocessor communication. 
Since typical FEM matrices are unstructured, algorithms such as the Recursive Spectral Bisection 
(RSB) have been found very effective in reducing inter-processor communication. However, the 
standard Gibbs-Pool-Stockmeyer profile reduction algorithm has been found even more effective 
in reducing the initial FE-ABC matrix (see Figure 3) to banded form as illustrated in Figure 10. 
This type of matrix reordering can deliver speed-ups as close to linear as possible. 



Figure 10: Reduced bandwidth of the FE-ABC system after application of the Gibbs-Pool-Stockmeyer 
profile reduction algorithm 
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5.3. Matrix Preconditioning 

Preconditioned iterative solvers are intended to improve the convergence rate of the algorithm. 
At times, preconditioners are necessary as may be the case with some dielect ricallv loaded struc- 
tures. However, for relatively small systems t less than 100.000 unknowns I it has been found that 
diagonal preconditioning is typically most effective and should always be applied. This precondi- 
tioning amounts to normalizing each row by the largest element, but even this simple operation 
can lead to substantial convergence speed-ups. Block and incomplete LF preconditioners are 
more effective in improving the convergence of the solver but are more costly to implement and 
one must judge on the overall CPF requirements rather than on the improved convergence alone. 
For example, the incomplete LI preconditioner given in [6Sj reduced the iterations to 1 /3 of those 
needed with diagonal preconditioning. However, each iteration was 3 times more expensive due 
to the triangular solver bottleneck. 


6. Reduced Order Modeling (ROM) of Frequency 
Responses 

Reduced Order Modeling (ROM) methods such as the Asymptotic Waveform Evaluation (AWE) 
have been successfully applied in V LSI and circuit analysis to approximate the transfer function 
associated with a given set of ports/variables in circuit networks [69, 70, 71]. The basic idea 
of the method is to develop an approximate transfer function of a given linear system from a 
limited set of spectral solutions. Typically, a Pade expansion of the transfer function is postulated 
whose coefficients are determined by matching the Pade representation to the available spectra] 
solutions of the complete system. 

In the context of finite element systems, ROM can be employed to predict a frequency response 
of the antenna input impedance or the scattering cross section of a given structure from a few 
data points. That is, once a few frequency points have been obtained by solving the entire finite 
element system of equations, these solutions along with the associated matrices can be re-used 
to extrapolate a broadband response without a need to resolve the system at other frequency 
points. In this section we present the theoretical basis of ROM and demonstrate its validity for 
full wave simulations using the finite element method as the computational engine. 

In addition to using ROM for antenna impedance and radar cross section prediction as a 
function of frequency, the method can also be used to fill-in angular pattern data points, thus 
eliminating a need to recompute the entire solution at small angular intervals. Since typical 
partial differential equation (PDE) systems involve thousands of unknowns, ROM can indeed 
lead to dramatic reductions of CPU requirements in generating a response of antenna or scatterer 
without a need to resolve the system for the fields in the entire computational grid. However, it 
should be noted that the FEM matrix for the reference frequency points must be stored (in core 
or out of core) with the current development of ROM for frequency domain analysis and thus 
some trade-off between CPU and memory requirements is unavoidable. Nevertheless, in view of 
the large CPU saving afforded by ROM, this appears to be a very small price to pay 
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6.1. Theory of Reduced Order Modeling 

6.1.1. FEM Syst( w 

When the functional (2) is discretized in connection with absorbing, boundarv conditions or 
artificial absorbers for mesh truncation, the resulting system can be decomposed into the form 



where A, denote the usual square (sparse) matrices and k = 2~/A = *c/c is the wavenumber of 
the medium. As usual. {/} is a column matrix describing the specific excitation. 

Clearly (24) can be solved using direct or iterative methods for a given value of the wavenum- 
ber as described earlier. Even though A, is sparse, the solution of the svstem (24) is computa- 
tionally intensive and must be repeated for each k to obtain a frequency response. Also, certain 
analyses and designs may require both temporal and frequency responses placing additional 
computational burdens and a repeated solution of (24) is not an efficient approach in generating 
these responses. An application of ROM to achieve an approximation to these responses is an 
attractive alternative. For these problems, the excitation column {/} is a linear function of the 
wavenumber and can therefore be stated as 

{/} = M/i} (25) 

with {/i } being independent of frequency. This observation will be specifically used in our 
subsequent presentation. 

6.1.2, Reduced Order Modeling 

To describe the basic idea of ROM in conjunction with the FEM, we begin by first expanding 
the solution {A"} in a Taylor series about Ar a , the wavenumber at which the system solution is 
available. We have 

{X} = {X a } + (k- k a ) {A,} + (k- k a ) 2 {X 2 } + ... 

+(* - *.)' {X,} + O {(k — fc a )' +1 } (26) 

where {A" a } is the solution of (24) corresponding to the wavenumber k a . By introducing this 
expansion into (24) and equating equal powers of k in conjunction with (25), after some manip- 
ulations, we find that 

{*.} = ULoM/i} 

{ X i} = A^ 1 [{f 1 }-A 1 {X a }-2k a A 2 {X a }} 

{X 2 } = -A- 1 [A 1 {A 1 } + A 2 ({A a } + 2fco{A 1 })] (27) 


{A,} = — Aq 1 [Ai + A 2 ({A/_ 2 } + 2k 0 {AMi})] 

with 

Aq = Aq T koAi + &qA 2 (28) 
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Expressions (27) are referred to as the system moments whereas 1 2 s ) is the system al the 
prescribed wavenumber ( Ay ). Although an explicit inversion of A" 1 may be needed as indicated in 
(27). this inversion is used repeatedly and can thus be stored out of-core for the implementation 
of ROM. Also, given that for input impedance computat ions we are typically interested in the 
field value at one location of the computational domain, only a single entrv of ]A',(A-)} need 
be considered, say (the pth entry) A 'f(k). The above moments can then be reduced to scalar 
form and the expansions (27) become a scalar representation of A'/’ (A) about the corresponding 
solution at A a . To yield a more convergent expression, we can instead revert to a Fade expansion 
which is a conventional rational function. 

For transient analysis the Pade expansion can be cast by partial fraction decomposition 
[71, 74] into 


-'?«••> = A To + Y, 


A* A*q A-, 


(29) 


where A', 0 is the limiting value as k tends to infinity. This is a qth order representation suitable 
for time/frequency domain transformation. As can be realized, the residues and poles (r, and 
k a + fa) in (29) correspond to those of the original physical system and play important roles 
in the accuracy of the approximation. As can be expected a higher order expansion with more 
zeros and poles can provide an improved approximation. The accuracy of ROM relies on the 
prediction of the dominant residues and poles located closest to Ay in a complex plane. Its key 
advantage is that for many practical electromagnetic problems only a few poles and zeros are 
needed for a sufficiently accurate representation. 

For a hybrid finite element - boundary integral system, the implementation of ROM is more 
involved because the fully populated boundary integral sub-matrix of the system has a more 
complex dependence on frequency. In this case we may instead approximate the full sub-matrix 
with a spectral expansion of the exponential boundary integral kernel to facilitate the extraction 
of the system moments. This approach does increase the complexity in implementing ROM. 
However, ROM still remains far more efficient in terms of CPU requirements when compared to 
the conventional approach of repeating the system solution at each frequency. 

As an application of ROM to a full wave electromagnetic simulation, we consider the evalu- 
ation of the input impedance for a microstrip stub shielded in a metallic rectangular cavity as 
shown in figure 11. The stub’s input impedance is a strong function of frequency from 1-3.2 GHz 
and this example is therefore a good demonstration of ROM’s capability. The shielded cavity is 
2.38cm x 6.00cm x 1.06cm in size and the microstrip stub resides on a 0.35cm thick substrate 
having a dielectric constant of 3.2. The stub is 0.79cm wide and A/2 long at 1.785 GHz and we 
note that the back wall of the cavity is terminated by a metal-backed artificial absorber having 
relative constants of e r = (3.2, -3.2) and /i r = (1.0, -1.0). 

As a reference solution, the frequency response of the shielded stub was first computed from 
1 to 3.2 GHz at 40 MHz intervals (50 points) using a full wave finite element solution. To 
demonstrate the efficacy and accuracy of ROM we chose a single input impedance solution at 
1.78 GHz. in conjunction with the 8th order ROM in (29) to approximate the system response. 
As seen in Figure 12, the 8th order ROM representation recovers the reference solution over the 
entire frequency band for both the real and reactive parts of the impedance. 

We conclude that the ROM technique is an extremely useful addition to electromagnetic 
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simulation codes and packages for computing wideband frequency response 
cross sect ionfRCS J pattern signature, etc. using only a feu samples of tin- 


's. large sets of radar 
syst cm solut ion. 



Figure 11: Illustration of the shielded microstrip stub excited with a current probe. 



Figure 12: Real and imaginary parts of input impedance computations based upon the 8th order ROM 
implementations using a single point expansion at 1.78 GHz. Solid lines: exact reference data; Dashed 
lines: 8th order ROM results. 


7. Additional Applications 

We choose two more examples to demonstrate the capability of the hybrid finite element methods. 

Scattering by a Large Cone-Sphere: A cone-sphere is basically a hemisphere attached to 
a cone. This is a difficult geometry to mesh since a surface singularity exists at the tip of the 
cone. The singularity can be removed in two ways: i) by creating a small region near the tip 
and detaching it from the surface or ii) by chopping off a small part near the tip of the cone. 
The second option inevitably leads to small inaccuracies for backscatter from the conical tip; 
however, we chose this option since the conical angle in our tested geometry was extremely small 
(around 7°) and the mesh generator failed to mesh the first case on numerous occasions. In 
Figure 13, we plot the backscatter patterns of a 4.5A long cone-sphere having a radius of 0.5A for 
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00 polarization. The mesh truncation surface is a rectangular Rox placed 0.1 A from the surface 
of the cone-sphere. As seen, the far-field results compare ext remely well with computations from 
a body of revolution code [75]. 
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Figure 13'. Backscatter pattern of a perfectly conducting conesphere for and 60 polarizations. 
Black dots indicate computed values using the FE-ABC code (referred to as FEMATS) and the solid 
line represents data from a body of revolution code [75]. Mesh termination surface is a rectangular box. 

Frequency Selective Surfaces (FSS): FSS structures [76] are arrays of tightly packed periodic 
elements which are typically sandwiched between dielectric layers. The periodic elements may 
be of printed form or slot configurations designed to resonate at specific frequencies. As such, 
they are penetrable around the element resonances and become completely reflecting at other 
frequencies. To meet bandwidth design specifications, stacked element arrays may be used in 
conjunction with dielectric layer loading. Here we consider the analysis of FSS structures via 
the FE-BI method. Because of the fine geometrical detail associated with the FSS surface, the 
finite element method has yet to be applied for the characterization of FSS structures, but use 
of prismatic elements makes this a much easier task. Of particular interest in FSS design is the 
determination of the transmission coefficient as a function of frequency, and since the array is 
periodic, it suffices to consider a single cell of the FSS. For computing the transmission coefficient 
T, the periodic cell is placed in a cavity as shown in Figure 14 and the structure is excited by a 
plane wave impinging at normal incidence. Assuming that near resonance the wave transmitted 
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Figure 14'. Upper figure: geometry of the multilayer frequency selective surface (FSS) used for modeling; 
lower figure: measured and calculated transmission coefficient through the FSS structure 


through the FSS screen will retain its TEM character, the transmission coefficient of the FSS 
panel can be approximated as 


1$ = 10 log 

where a is the reflection coefficient of the absorber placed at the bottom of the cavity and should 
be kept small (< 0.1) to suppress higher order interactions. By adding the next higher order 
interaction, a more accurate expression for the transmission coefficient is 

TdB ~ Tjg + 10 log [l-a(l-T<°>)]. 

The above FSS modeling approach was applied for a characterization of multi-layered slot 
FSS structures. The geometry of the multilayer radome is given in Fig. 14. The total thickness 
of the FSS was 6.3072cm and is comprised of two slot arrays (of the same geometry) sandwiched 
within the dielectric layers. For modeling purpose, a 1.54cm thick absorber is placed below the 
FSS as shown in Fig 14. It is seen that the results generated by the FE-BI method are in good 
agreement with the measurements [77]. 


a 
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8. Conclusion 

WV- reviewed hybrid finite element methods as applied to electromagnetic scattering and radia- 
tion problems. Much of the emphasis dealt with the various mesh truncations schemes and we 
presented an up-to-date account of these schemes. The usual finite element -boundary integral 
method was presented and new developments for reducing the CRT requirements of this tech- 
nique using the fast integral methods were discussed. Antenna feed modeling in the context 
of the finite element method had not been discussed before and for the first time we presented 
an overview of the modeling approaches for the most popular antenna feeds, including aperture 
coupled feeds. Parallelization will continue to play an increasingly greater role and a section 
was included discussing our experiences for better implementation of finite element codes on 
distributed and vector architectures. A number of examples illustrating the successful applica- 
tion of the finite element method were included throughout the paper and these were intended 
to demonstrate the method's geometrical adaptability and inherent capability to treat highly 
heterogeneous structures. 

As can be expected, issues relating to mesh truncation, mixing of elements [78], domain 
decomposition [79, 80], robustness, adaptive refinement[81], accuracy, error control, feed mod- 
eling and parallelization for large scale simulations will continue to dominate future research 
and developments relating to partial differential equation methods. Reduced order modeling 
techniques such as the AWE method are also very promising for reducing the computational re- 
quirements in generating broadband responses. Further development of AWE is certainlv needed 
for its application in connection with hybrid finite element systems. 

An apparent advantage of the finite element method is its potential hybridization with all 
other frequency domain methods. Future applications of the finite element method are likely 
to make greater use of hybridization techniques aimed at increasing the method’s accuracy and 
efficiency while retaining its inherent geometrical adaptability and ease in handling materials. 
Reduced order modeling techniques such as the AWE method is another promising approach 
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